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SUMMARY 


Two  mathematical  models  that  simulate  sea  ice  motions  from  given 
wind  and  current  inputs  are  described.  These  models  are  evaluated  as 
candidates  for  predicting  and  simulating  oil  spill  trajections  in  ice- 
covered  waters.  Simulated  displacements  of  ice  are  compared  with  observed 
displacements  to  evaluate  the  performance  of  these  models  during  the 
AIDJEX  main  field  experiment  in  1975-76.  During  this  year-long  experi¬ 
ment,  accurate  winds  and  currents  were  available  as  input  to  the  ice 
motion  models  and  a  set  of  data  buoys  and  camps  provided  motion  observa¬ 
tions  for  evaluating  model  performance.  In  addition  to  this  performance 
validation,  requirements  for  computer  resources,  manpower  and  input  data 
are  also  presented. 

A  free-drift  ice  motion  model  is  studied  as  the  simplest  possible 
model  for  simulating  daily  ice  motions.  In  this  model,  forces  due  to 
winds,  currents,  Coriolis  acceleration  and  sea  surface  tilt  are  balanced. 
The  daily  ice  motion  at  any  selected  site  is  obtained  from  the  local 
wind  and  current  during  that  day.  During  the  summer  of  1975,  say  May 
through  September,  simulated  daily  ice  motions  are  compared  with  observed 
motions.  The  mean  displacement  error  is  0.010  m/s  and  the  standard 
deviation  is  0.030  m/s.  These  errors  may  be  attributed  to  uncertainty 
in  winds  and  currents.  During  the  winter  of  1975-76  when  ice  stress  is 
important,  these  errors  increase.  Then,  a  mean  displacement  error  of 
0.030  m/s  and  standcrd  deviation  of  0.107  m/s  are  apparent.  Errors 
increase  roughly  linearly  during  the  fall,  say  October  through  December, 
and  decrease  rapidly  during  April  and  May  (of  1975,  the  previous  year). 

A  complete  ice  dynamics  model  that  describes  internal  ice  stress 
according  to  an  elastic  plastic  material  description  is  also  studied. 

This  complete  model  is  used  to  simulate  ice  motions  at  three  selected 
times  during  the  study  period  in  1975-76.  During  summertime  and  early 
fall  the  complete  ice  dynamics  model  is  slightly  more  accurate  than  the 
free-drift  model,  with  mean  displacement  errors  being  0.005  m/s  and 
standard  deviation  being  0.020  m/s.  Although  slightly  better,  these 
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errors  are  comparable  with  those  of  free  drift.  However,  dramatic 
differences  in  the  size  of  errors  appear  in  late  winter  (January  and 
February  1976).  During  this  time  the  complete  ice  dynamics  model  gives 
•mean  displacement  errors  of  0.005  m/s  and  a  standard  deviation  of  0.036 
m/s.  From  this  evidence  we  are  able  to  conclude  that  the  complete  ice 
dynamics  model  is  capable  of  describing  ice  motions  to  within  the  same 
accuracy  throughout  the  year. 

Both  mathematical  models  of  sea  ice  dynamics  require  winds  and 
currents  to  determine  the  ice  motions.  The  complete  ice  model  also 
requires  ice  conditions  and  far-field  boundary  motion  information.  The 
free-drift  computer  requirements  are  minimal,  while  the  complete  model 
requires  a  computer  with  roughly  100K  octal  words  of  central  memory. 
However,  the  input  of  wind  and  current  fields  for  both  models  probably 
requires  this  size  of  computing  facilities  and  graphic  output  of  the 
results  also  dictates  that  a  substantial  computer  be  used.  It  is  at 
present  felt  that  manpower  requirements  for  operational  forecasting  are 
similar  for  both  models. 

Both  the  free-drift  and  complete  ice  dynamics  models  are  useful 
tools  for  U.S.  Coast  Guard  requirements  of  predicting  motions  of  oiled 
ice  in  the  Beaufort  Sea.  During  wintertime  the  complete  model  must  be 
used  if  maximum  accuracy  is  to  be  obtained,  but  during  summertime  the 
free-drift  model  appears  to  be  nearly  as  accurate  as  the  complete  model. 
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Chapter  1 
INTRODUCTION 


1.1  OBJECTIVE 


Recent  lease  sales  (December  1979)  of  offshore  tracts  in  the  Beaufort 
Sea  portion  of  the  Arctic  Ocean  represent  the  beginnings  of  the  development 
cf  oil  fields  in  this  marine  environment.  Exploratory  drilling  operations 
in  both  the  U.S.  and  Canada  must  be  expected  to  increase  and  be  followed 
by  production  operations  as  the  oil  fields  are  developed.  The  thick  sea 
ice  cover  on  the  Arctic  Ocean  is  a  formidable  factor  in  these  operations 
Mobile  pack  ice  is  present  year-round  in  the  central  Arctic.  It  is 
composed  of  many  thick,  multiyear  ice  floes  that  are  typically  a  few 
kilometers  in  si2e  and  are  separated  by  thin  ice  more  recently  grown  on 
leads  formed  as  the  floes  fracture  and  divide.  The  large-scale  average 
thickness  of  the  pack  ice  is  roughly  3  m.  As  an  opposite  extreme, 
during  the  summer  when  air  temperatures  rise  above  freezing,  melting  and 
offshore  winds  combine  to  form  an  approximately  300-km-wide  swath  of 
open  water  along  the  Alaskan  North  Slope  and  the  Mackenzie  Delta.  Thus, 
summertime  drilling  could  take  place  (U.S.  regulations  do  not  allow 
summertime  drilling  in  the  Prudhoe  Bay  lease  area,  but  it  is  the  only 
time  drilling  is  allowed  in  the  Mackenzie  Bay  area  of  Canada)  using 
techniques  developed  for  shallow,  midlatitude  waters.  After  freeze-up, 
the  coastal  or  continental  region  is  ice  covered  except  for  leads  and 
polynyas  generated  by  ice  deformations.  The  ice  nearest  shore  is  typically 
motionless  because  it  is  often  grounded,  is  protected  by  shoreline 
irregularities  and  is  strong  enough  to  resist  deformation  by  winds  and 
currents.  Between  the  mobile  pack  ice  and  the  land-fast  ice  is  the 
highly  deformed  stamukhi  zone,  an  area  in  which  shear  and  pressure 
ridges  are  found  nearly  everywhere. 

This  quick  sketch  of  the  ice  conditions  that  appear  in  Arctic 
waters  is  included  here  to  help  indicate  the  range  of  applicability  of 
the  models  presented  in  this  report.  Variations  in  the  location  and 
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extent  of  each  ice  zone  are  significant  from  year  to  year  and  during  the 
course  of  each  year.  The  large-scale  ice  drift  models  are  intended  to 
describe  the  behavior  of  pack  ice,  not  the  smaller-scale  fast  ice. 

Thus,  some  care  must  be  exercised  to  ensure  that  the  models  are  applied 
only  under  conditions  which  adequately  simulate  the  appropriate  physical 
processes . 

This  study  is  part  of  a  U.S.  Coast  Guard  program  designed  to  provide 
methods  to  predict  the  movement  of  oil  spills  in  Arctic  and  sub-Arctic 
coastal  waters.  These  motion  predictions  are  necessary  to  allow  oil 
spill  cleanup  strategies  to  be  selected  and  plans  to  be  developed  for 
handling  individual  oil  spills.  Defining  the  kinds  of  spills  that  can 
occur  is  not  part  of  this  study,  but  it  is  useful  to  note  that  an  oil 
spill  could  range  from  a  small  release  of  fixed  volume  to  an  oil  well 
blowout  that  continues  unchecked  for  an  entire  year.  Thus,  a  method  is 
needed  for  tracking  oil  spills  over  periods  ranging  from  a  few  days  to 
many  months. 

A  great  amount  of  research  has  been  performed  to  determine  the  fate 
and  behavior  of  oil  spilled  in  and  under  a  sea  ice  cover.  A  synthesis 
of  these  results  has  been  prepared  as  part  of  the  Outer  Continental 
Shelf  Environmental  Assessment  Program  (Thomas,  19t0) .  The  most  im¬ 
portant  result  of  that  work  is  that  oil  will  be  transported  primarily  by 
the  large-scale  motions  of  the  pack  ice.  Therefore,  tracking  of  an  oil 
spill  becomes  a  problem  of  tracking  the  ice  cover  itself.  The  work 
presented  in  this  study  is  focused  on  evaluating  the  performance  of 
models  that  simulate  pack  ice  motions  to  determine  probable  ice  trajectories. 

1 . 2  APPROACH 

Two  different  mathematical  models  of  sea  ice  dynamics  are  evaluated 
to  learn  how  well  each  simulates  observed  motions  under  a  full  range  of 

'The  Outer  Continental  Shelf  Environmental  Assessment  Program  (OCSEAP) 
is  funded  by  the  Bureau  of  Land  Management  through  interagency  agree¬ 
ment  with  the  National  Oceanic  and  Atmospheric  Administration  to  respond 
to  the  needs  of  petroleum  development  off  the  Alaskan  Continental  Shelf. 


2 


ice  conditions.  The  first  is  a  free-drift  model  in  which  air  and 
water  drag,  Coriolis  force  and  sea  surface  tilt  are  the  forces  acting 
on  the  ice.  This  model  neglects  the  ice  stress  that  arises  as  a  larger- 
scale  effect  of  forces  transmitted  between  individual  ice  floes.  It  is 
the  simplest  deterministic  formulation  that  can  be  derived  based  on  the 
driving  forces  applied  by  the  atmosphere  and  the  ocean.  The  second  is 
a  complete  ice  dynamics  model  that  describes  the  internal  ice  stress  as 
a  hardening/sof tening,  elastic-plastic  material  in  which  strength  is  a 
function  of  the  instantaneous  ice  thickness  distribution.  This  model, 

X 

developed  during  AIDJEX  ,  is  the  most  sophisticated  mathematical  repre¬ 
sentation  of  the  physical  processes  that  control  sea  ice  dynamics  behavior 
that  has  been  developed.  The  major  features  of  the  AIDJEX  model  are 
presented  by  Pritchard  (1980) ;  they  are  modified  extensively  and  reformu¬ 
lated  in  his  later  work  (Pritchard,  1981). 

The  performance  of  each  model  is  evaluated  by  determining  the 
difference  between  observed  and  simulated  daily  displacements  (average 
velocities  during  a  day)  using  data  from  AIDJEX  buoys  and  manned  camps. 
Several  time  periods  during  the  1975-76  AIDJEX  main  experiment  are 
selected:  13-28  April  1975,  17-25  May  1975  and  27  January  through 
12  February  1976.  These  are  the  only  periods  for  which  AIDJEX  model 
simulations  are  available.  In  addition,  free-drift  simulations  are 
evaluated  from  June  through  December  1975.  The  scope  of  the  work  presented 
in  this  report  is  limited  to  evaluating  the  performance  and  estimating  the 
resources  required  to  use  each  model  for  oil  spill  tracking.  However,  our 
primary  focus  in  this  report  is  to  evaluate  the  free-drift  model  and  to 
learn  how  well  it  performs  in  comparison  to  the  more  complicated  model. 

The  free-drift  model  is  described  completely.  In  this  discussion, 
all  forces  acting  on  the  ice  (except  internal  ice  stress)  are  presented. 

The  planetary  boundary  layer  model  and  oceanic  mixed  layer  models  show 


‘The  Arctic  Ice  Dynamics  Joint  Experiment  (AIDJEX),  1970-77,  was  a  multi- 
nationally  coordinated  theoretical  and  field  program  aimed  at  gaining  a 
better  fundamental  understanding  of  the  physical  behavior  of  sea  ice. 
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how  traction  is  transmitted  from  winds  and  currents.  The  minimum  time 
resolution  of  1  day  is  identified  and  the  ice  response  is  explained. 
The  complete  ice  dynamics  model,  which  is  a  reformulation  of  the  early 
AIDJEX  model,  is  also  presented  in  detail.  The  input  data  and  the 
computing  resources  required  to  use  each  model  are  listed. 
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Chapter  2 
FREE-DRIFT  MODEL 


e 


2.1  GENERAL  DESCRIPTION 

In  1979  the  Outer  Continental  Shelf  Environmental  Assessment  Program 
(OCSEAP)  funded  a  study  by  Flow  Research  Company  of  ice  motions  in  the 
Beaufort  and  Chukchi  Seas  (Thomas  and  Pritchard,  1979).  The  free-drift 
model  as  taken  from  that  report  is  described  in  this  section.  The 
effect  of  model  parameters  on  free-drift  is  studied  by  nondimensionalizing 
the  force  balance  equations  driving  the  model. 

An  important  consideration  in  applying  the  free-drift  model  is  time 
scale.  Inertial  effects  are  ignored  in  free-drift  so  that  the  model  has 
a  resolution  of  1  day.  This  is  short  enough  to  account  for  synoptic 
variations  in  weather  systems  while  still  being  long  enough  to  average 
out  high  frequency  motions.  If  an  ice  trajectory  for  a  month  is  desired, 
it  is  tempting  to  apply  the  wind  stress  averaged  over  a  month  and  make  a 
single  calculation.  This  is  not  acceptable  because  nonlinearities  in 
the  model  require  daily  calculations  of  drift. 

If  a  free-drift  model  is  to  be  implemented,  it  is  important  to  know 
the  quality  of  the  data  required  to  drive  it.  Free-drift  calculations 
require  both  air  and  water  stress  data  for  daily  intervals  to  be  input. 
Both  of  these  parameters  are  difficult  to  estimate  accurately.  Estimates 
of  the  air  stress  require  knowledge  of  the  wind  field,  which  is  generally 
calculated  from  the  distribution  of  atmospheric  pressure.  In  the  Arctic 
Ocean,  this  distribution  is  not  well  known.  The  water  stress  is  due  to 
the  relative  motion  of  the  ice  with  respect  to  the  water.  Estimates  of 
water  stress  thus  require  knowledge  of  geostrophic  currents.  These 
currents  vary  from  0.01  to  0.03  m/s  in  the  central  Beaufort  and  are  not 
well  known  even  for  1975  when  many  measurements  were  made.  In  practice, 
and  for  all  the  free-drift  calculations  presented  here,  it  is  necessary 
to  assume  a  zero  geostrophic  current  and  accept  this  as  a  major  source 
of  error. 
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2.2  FREE-DRIFT  EQUATIONS  OF  MOTION 


During  free-drift,  motion  of  the  ice  cover  may  be  determined  by 

coosidering  momentum  balance  locally.  The  forces  acting  on  the  ice  are 

air  stress  T  ,  water  stress  x  ,  Coriolis  force  -mf  kxv  and  sea 
~a  ~w  ~  ~ 

surface  tilt  -nvjVH.  Momentum  changes  then  occur  as 

mv  =  T  +  T  -  mf  kXv  -  rmcVH  (2.1) 

~  ~a  ~w  ~  ~ 

where 

m  is  mass  per  unit  area  of  ice, 

v  is  velocity  in  the  horizontal  plane, 

v  is  material  time  rate  of  change  of  v, 

is  traction  exerted  by  the  atmosphere  on  the  upper  surface  of  ice, 
x  is  traction  exerted  by  the  ocean  on  the  lower  surface  of  ice, 
f  is  the  Coriolis  parameter, 

k  is  the  unit  vector  in  the  vertical  direction, 
g  is  acceleration  due  to  gravity,  and 
H  is  the  height  of  the  sea  surface. 

The  air  stress  is  determined  from  the  geostrophic  wind  in  the 
atmosphere  as 


T 

~a 


P  c 
a  a 


|U|  ?aU 


(2.2) 


where 


'  cos  a  -sin  a 


B 

~a 


ism  a  cos  a 


(2.3) 


and  p  is  air  density  and  c  is  the  drag  coefficient.  The  air  stress 

3  3. 

is  turned  at  an  angle  a  to  the  applied  geostrophic  wind  U  .  We  have 
neglected  the  ice  drift  since  it  is  two  orders  of  magnitude  smaller  than 
the  geostrophic  wind.  The  geostrophic  wind  is  determined  from  the 
surface  barometric  pressure  field  P  as 

U  -  -  tO(VP  .  (2.4) 
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The  water  stress  satisfies  a  quadratic  drag  law  similar  to  the 


atmosphere: 


»  T=pc|v-v|B(v-v)  (2.5) 

~w  w  w  ~g 1  ~  ~  ~g 

where  p  is  water  density,  c  is  the  drag  coefficient, 
w  *  w 

ocean  geostrophic  current  and 

sin  (8+tt)  cos  (8+tt)/ 

where  8  is  the  angle  at  which  traction  is  applied  by  the  current. 

The  sea  surface  tilt  defines  the  ocean  geostrophic  current  in  the  form 

mcVH  =  -  mf  kXv  .  (2.7) 

- g 

The  ice  velocity  may  be  determined  at  each  point  as  a  function  of 
time  whenever  the  barometric  pressure  field  history  is  prescribed.  Many 
authors  have  used  the  free  drift  relationships  to  analyze  ice  motions. 
Recently,  McPhee  (1980)  has  used  this  model  to  simulate  motion  and 
deformation  of  the  AIDJEX  manned  array  during  the  summer  of  1975. 

The  results  sought  have  a  time  resolution  on  the  order  of  1  day. 

For  this  case,  inertia  is  negligible.  Therefore,  the  analysis  is 

performed  for  steady-state  conditions.  When  shorter  time-scale  results 

are  desired,  the  ice  inertia  may  be  included;  but  in  this  case  the 

inertia  of  the  water  column  is  also  significant  (McPhee,  1978).  Thus, 

the  inertial  contribution  is  not  completely  formulated.  For  both 

reasons  we  choose  to  neglect  inertia  completely  and  use  the  steady-state 

equations  to  calculate  daily  average  velocities.  In  this  case.  Equation  (2.1) 

becomes 


v  is  the 
~g 


(2.6) 


Ta  +  Tw  -  mf  kXv  -  mgVH  =  0  (2.8) 

If  the  water  stress  (2.5)  and  sea  surface  tilt  terms  (2.7)  are  substituted 
into  Equation  (2.8),  then 


7 


(2.9) 


t  +  pc  |v-v|B(v-v)-mfkx(v-v)=0 
~a  Kw  w  1 ~  ~g'  ~  ~  ~g  -  ~  -g 


BuC  this  force  balance  relates  air  stress  T  ,  geostrophic  ocean  currents 

■—  3 

v.  and  ice  drift  in  a  special  way,  namely,  the  ice  velocity  and  current 
”  S 

always  appear  as  the  difference  v  -  v  .  Thus,  if  we  introduce 


G  = 


(2.10) 


then  force  balance  becomes 

T  +  p  c  I G I  BG  -  mf  kxG  =  0 

This  equation  may  be  rewritten  in  terms  of  water  stress  if  desired  as 


(2.11) 


T  +  T  -  of  icxG  =  0 
-a  ~w  ~  ~ 


(2.12) 


where 


T  *  P  C  ;G(  BG  . 
~w  w  w  ■ ~  — 


As  a  result  of  introducing  G  ,  the  ice  velocity  relative  to  deep 
ocean  currents,  it  is  seen  that  it  is  this  quantity  that  air  stress 
affects.  Therefore,  contributions  to  ice  drift  due  to  air  stress  and  to 
ocean  current  may  be  calculated  separately.  The  air  stress  effect  is 
found  from  Equation  (2.11)  and  the  ocean  current  is  added  directly  in 
the  form 


v  =  G  +  v 
-g 


(2.13) 


2.3  MODEL  RESPONSE 

In  order  to  understand  model  response,  nondimens ional  variables  are 
introduced  that  show  explicitly  how  all  material  constants  affect  solutions 
and  simplify  the  presentation  and  understanding  of  the  results.  The 
nondimensional  variables  have  been  determined  by  Pritchard  and  Thomas  (1980) 
for  the  AIDJEX  ice  model. 
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We  introduce  nondimens ional  variables  for  velocity, 


P  c 

—  W  w 

V  =  - T—  V 

-  mf 


P  c 

—  w  w 

u  =  u 

-  mf 


P  c 

?  =  -w  w  r 
5  mf  5  ’ 


and  for  air  and  water  stress, 


—  PwCw 

1  =  - 2  I 

(mf  j 


and  for  pressure  and  position  we  introduce  arbitrary  scaling. 


(2.14) 


(2.15) 


p  -  p/ps 

x  =*  x/x^ 


We  introduce  the  nondimens ional  numbers 


N 


p  c  P„ 
w  w  R 


P  pa 


P  c 

N  = 
a  pc 
w  w 


and  the  steady  state  equation  of  motion  (2.11)  then  becomes 


X  +  I G I  BG  -  kXG  *  0 
~a  1  - 1  ~~  ~  - 

where  air  stress  satisfies 


X  =  N  |  U I  B  U  . 
-a  a  ~ 1  ~a~ 

The  geostrophic  wind  is 

U  -  -N  kxV  P 

p  - 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


(2.20) 
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The  relationship  between  ice  velocity  and  applied  air  stress  has 
been  described  in  a  convenient  form  by  Pritchard,  Coon  and  McPhee  (1977). 
If  we  modify  that  work  by  using  the  present  nondimensional  variables, 
then 


where 


sQG 


=  (  !  G  !  2  +  2  1 G  |  sin  £  +  1 


cos  o  -sin  $ 


.sin  6  cos  c 


(2.21) 


(2.22) 


where 


tan  5  *  tan  £  +  — — -  .  (2.23) 

|G(  cos  £ 

The  ice  drift  relative  to  geostrophic  ocean  currents  is  to  the  right  of 
the  air  stress  by  the  angle  i  .  The  forces  acting  on  the  ice  are  shown 
in  Figure  2.1.  The  ice  drift  appears  as  a  dashed  arrow. 


Figure  2.1.  Free-Drift  Force  Balance  and  Ice  Velocity 
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The  magnitude  of  air  stress  and  ice  drift  are  related  by 

ItJ  -  s  |G|  ,  (2.24) 

which  is  displayed  in  Figure  2.2.  We  have  chosen  to  invert  Equation  (2.24) 
for  graphical  representation  since  we  anticipate  T  to  be  given  as 

*-  3 

input  and  ice  drift  calculated  as  output.  The  angle  5  also  varies 
with  air  stress  magnitude  as  may  be  seen  by  solving  Equations  (2.23)  and 
(2.24)  simultaneously.  In  Figure  2.3  we  present  the  results.  It  is 
seen  that  the  turning  angle  depends  significantly  on  the  oceanic  boundary 
layer  Ekman  angle  8  .  This  is  not  true  of  ice  speed  which  is  insensitive 
to  changes  in  8  • 

The  ice  speed  is  a  nonlinear  function  of  air  stress  magnitude, 
and  the  turning  angle  depends  on  air  stress  magnitude,  too.  Thus,  the 
response  is  nonlinear  and  therefore  the  model  cannot  be  used  directly  on 
monthly  mean  air  stress. 

To  circumvent  this  problem,  we  look  at  the  relationship  between  air 
stress  magnitude  and  ice  drift  speed  which  is  nearly  parabolic.  By 
introducing  the  quadratic  atmospheric  boundary  layer  model,  we  seek  a 
linear  relationship  between  geostrophic  wind  and  ice  drift.  We  find 

N  | U |  B  U  =  sQG  .  (2.25) 

The  wind  and  ice  speed  are  related  by 

Na  |U|2  -  s  |G|  .  (2.26) 

This  relationship  is  presented  in  Figure  2.4.  It  is  seen  that  the  ice 
speed  is  nearly  a  linear  function  of  the  wind  speed.  This  observation 
was  noted  by  both  M.  G.  McPhee  (1980)  and  A.  S.  Thorndike  (personal 
communications)  and  led  them  to  prefer  using  geostrophic  wind  rather 
than  air  stress  in  their  free-drift  analyses. 
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Figure  2  3  Relative  Ice  Velocity  Orientation  <)  as  a  Function  of  Nondimensional  Air  Stress 
Magnitude  ItJ  for  a  Range  of  Oceanic  Boundary  Layer  Turning  Angles  fi 


The  angle  between  geostrophic  wind  and  ice  drift  is  controlled  by 
T 

B  Q  ,  where  the  superscript  T  represents  the  transpose  of  the  matrix 

B  and  since 
-a 

N  |U|  U  =  sBTQG  ,  (2.27) 

we  see  that  the  ice  drift  is  oriented  at  an  angle  5  -  a  to  the  right 
of  the  geostrophic  wind.  This  relationship  is  shown  in  Figure  2.5.  It 
is  seen  that  the  turning  angle  depends  strongly  on  wind  speed.  Therefore, 
the  relationship  is  not  linear  and  the  monthly  average  geostrophic  wind 
cannoc  be  used  directly  to  obtain  the  monthly  ice  motion.  In  practice, 
however,  it  is  possible  that  a  constant  angle  (say  S  =  45°)  would  provide 
a  useful  approximation  to  longer-term  drift. 

Due  to  the  nonlinearity  in  these  relationships,  the  decision  was 
made  to  process  historical  winds  on  a  daily  basis  and  sum  up  the  daily 
ice  displacements  to  obtain  monthly  values.  The  nondimens ional  formulation 
has  simplified  the  task,  of  determining  the  effect  of  material  parameters, 
however.  To  estimate  these  effects  better,  consider  the  relationship 
between  ice  speed  and  geostrophic  wind.  The  ice  speed  is  nearly  a 
linear  function  of  geostrophic  wind  speed  at  higher  values.  If  Equation  (2.26) 
is  rewritten  as 


2  /  2  . 
y  =  (x  +  ax  + 


x 


(2.28) 


where 

y  =  N  ^  |  U  j  ,  x  =  j  G  |  and  a  =  2  sin  3 

are  introduced  for  notational  convenience,  then  an  expansion  may  be  made 
for  large  values  of  x  : 


y 


(1  +  ax 


-2  ^ 
x  )  *  x 
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Ice  Turning  Angle,  d(degrees) 


Figure  2.5.  Relative  Ice  Velocity  Orientation  d  as  a  Function  of  Nondimensional 
Geostrophic  Wind  Speed  Ng2  ju|  for  a  Range  of  Oceanic  Boundary 
Layer  Turning  Angles  ft 


and  then 


y 


(1  +  ^  ax  *  +  . . . )x 


Finally, 


y 


=  J  +  x  +  0 (x-1) 


(2.29) 


Returning  to  the  original  notation  and  inverting  results  in 


|G| 


j  sin  S 


(2.30) 


which  provides  the  coefficients  of  the  linear  relationship.  If  the 
nondimensional  solutions  given  by  Equations  (2.23)  and  (2.30)  are  re¬ 
written  in  terms  of  physical  variables,  then 


and 


mf 


2p  c 
w  w 


sin  8 


(2.31) 


tan  6  =  tan  3  +  - - -  .  (2.32) 

p  c  G  cos  3 
w  w  ~ 1 

It  is  seen  that  the  ice  speed  is  insensitive  to  changes  in  mass  m 

because  it  occurs  only  in  the  constant  term,  which  is  small.  But  the 

ratio  of  drag  coefficients  in  the  atmosphere  and  ocean,  N  ,  affects 

a 

ice  speed  as  a  square  root  function.  Therefore,  all  parameters  have  a 

modest  affect  on  ice  speed.  However,  in  the  range  of  parameters  anticipated, 

all  can  cause  large  changes  in  the  turning  angle.  For  example,  increasing 

8  from  20°  to  30°  increases  the  turning  angle  by  about  10°.  This  is 

comparable  to  increasing  the  water  drag  or  decreasing  the  mass  density 

m  by  about  10  percent  at  higher  speeds.  The  changes  are  smaller  when 

winds  are  lower.  The  ratio  of  drag  coefficients,  N  ,  affects  the 

a 

turning  angle  indirectly  through  the  ice  speed.  If  a  10  percent  increase 
in  causes  a  5  percent  increase  in  [ G |  ,  this  in  turn  will  cause 

about  a  5°  decrease  in  the  turning  angle. 


17 


Some  consideration  must  be  given  to  the  special  conditions  that 
exist  near  the  edge  of  the  pack  ice  during  summer  when  there  is  appreciable 
open  water.  In  this  case,  the  ice  concentration  is  often  low.  The 
model  has  not  been  tested  thoroughly  under  these  conditions;  however,  the 
following  argument  should  eliminate  serious  concern.  In  the  free-drift 
model  (Equation  2.11),  the  areal  mass  density  is  scaled  down  in  proportion 
to  the  ice  concentration.  Similarly,  the  traction  exerted  by  the  atmosphere 
and  the  ocean  on  the  ice  in  Equations  (2.2)  and  (2.11)  should  be  proportional 
to  the  relative  area  covered  by  ice  -  but  this  is  just  the  concentration. 
Therefore,  each  term  in  the  free-drift  momentum  balance  is  changed  by 
the  same  amount,  and  the  relationship  between  ice  motion  and  geostrophic 
wind  is  unchanged. 

2 . 4  PARAMETERS 


The  main  driving  force  of  the  free-drift  model  is  the  air  stress. 
The  air  stress  is  computed  from  the  geostrophic  wind  field  as 


T 

~a 


B  = 
~a 


paca  'H>  5aH 


cos  a  -sin  a' 


sin  a  cos  a , 


(2.2  bis) 


(2.3  bis) 


where  c  is  a  drag  coefficient,  p  is  the  air  density,  U  is  the 
geostrophic  wind  and  a  is  the  angle  of  the  air  stress  from  the  geostrophic 
wind.  The  geostrophic  winds  are  computed  from  the  surface  pressure 
field  as 


U  - — r  !tX7P 

Paf  . 

where  f  is  the  Coriolis  parameter. 


(2.4  bis) 


18 


to  be 


Leavitt  et  ai.  (1978a, b)  have  estimated  c^q  ,  the  10-m  drag 
coefficient,  to  be  0.0025  using  AIDJEX  data.  They  believe  c^q 
accurate  to  +20  percent.  Albright  (1980)  has  computed  the  ratio  U^q/U 
to.be  0.55  in  the  summer,  using  the  AIDJEX  data  set.  Computing  the  drag 
coefficient  as 


and  using  the  above  values  of  U^q/U  and  c^q  gives  =  0.0008  for 
the  summer.  Albright  (1980)  also  gives  the  summer  value  of  a  =  24°  with 
a  standard  deviation  of  5°. 

The  air  density,  p^  ,  is  computed  as  the  monthly  average  of  all 
the  densities  computed  from  AIDJEX  temperature  data.  The  "within  month" 
standard  deviation  of  is  about  2  percent  while  the  range  of  variation 

from  month  to  month  is  about  12  percent. 

Reliable  data  on  the  mean  thickness  of  the  ice  for  1975  are  not 
available.  We  used  the  computed  mean  thickness  presented  in  Thorndike 
et  al.  (1975),  where  an  initial  ice  thickness  distribution  as  measured  by 
under-ice  profile  is  followed  through  2  years  of  known  deformation  and 
climatological  thermodynamics.  Table  2. 1  gives  the  air  density  and  ice 
thickness  used  in  free-drift  simulations.  One  check  point  of  these  data 
appears  during  April  1976  when  an  under-ice  profile  was  obtained  by  upward- 
looking,  submarine-mounted  sonar.  Analysis  of  these  data  by  Wadhams  and 
Horne  (1978)  shows  a  mean  draft  of  3.7  m  (average  of  two  tracks  of  200  km 
length  crossing  under  camp  Caribou).  Converting  mean  draft  to  mean  ice 
thickness  is  accomplished  by  multiplying  by  the  ratio  of  ocean-water 
density  to  ice  density  p/p^  “1.1.  The  mean  ice  thickness  during 
April  1976  is  observed  to  be  4.1  ra.  This  observed  value  is  about  5  percent 
larger  than  the  April  values  calculated  by  Thorndike  et  al.  (1975).  It 
is  likely  that  values  in  Table  2. 1  for  January  through  April  are  somewhat 
low  since  mean  thickness  should  increase  throughout  winter.  The  history 
shown  is  probably  dependent  on  initial  conditions  and  on  the  specific 
deformation  history  used  by  Thorndike  et  al.  (1975).  In  any  case,  the 
values  are  reasonable  and  an  uncertainty  of  5  percent  is  acceptable. 


19 


Table  2.1.  SEASONAL  VARIATIONS  IN  ICE  THICKNESS  AND  AIR  DENSITY. 


Month 

h  (m) 

Pa  0<-g/' 

January 

3.90 

1.47 

February 

3.90 

1.47 

March 

3.90 

1.45 

April 

3.90 

1.40 

May 

3.70 

1.36 

June 

3.55 

1.31 

July 

3.25 

1.29 

August 

2.95 

1.29 

September 

3.30 

1.33 

October 

3.70 

1.38 

November 

4.15 

1.44 

December 

4.05 

1.46 

The  Coriolis  parameter  is  computed  as 

f  =  2ft  sin  <J>  (2.  34) 

_5 

where  ft  ■  7.29  x  10  rad/s  and  $  =  latitude. 


The  water  stress  is  computed  as 


T 

~w 


pc  Iv-vlB(v-v) 
w  w  -g1  -  -  ~g 


(2.3  bis) 


(cos  (6+tt)  -sin  (6+tt) 

sin  (S+tt)  cos  (8+it) 


(2.6  bis) 


The  water  drag  coefficient  and  turning  angle  used  were  c  =  0.0055  and 

V 

6  *  23  ,  as  given  in  McPhee  (1980).  The  geostrophic  current  v 

— 8 

was  set  to  zero  as  previously  explained. 
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Chapter  3 

COMPLETE  ICE  DYNAMICS  MODEL 


3,1  GENERAL  DESCRIPTION 

The  sea  ice  dynamics  model  describes  the  large-scale  mechanical 
behavior  of  the  thin  ice  cover  floating  on  a  polar  ocean.  Thus,  we 
require  a  mathematical  description  of  forces  acting  on  a  large-scale 
typical  element  of  ice.  The  model  is  two-dimensional  in  the  sense  that 

all  motions  are  constrained  to  lie  in  the  surface  of  the  ocean.  But  the 

thickness  of  pack  ice  is  variable,  being  made  up  of  thick  old  floes, 
open  water  and  thin  young  ice  growing  on  leads.  The  amounts  of  open 
water  and  thin  ice  have  a  dominant  influence  both  on  the  rate  of  vertical 
heat  transfer  from  the  ocean  to  the  atmosphere  and  on  the  strength  or 
ability  of  the  ice  to  resist  wind  or  ocean  current  loads.  The  variations 
in  thickness  are  characterized  by  a  thickness  distribution  that  describes 
the  fraction  of  an  elemental  region  covered  by  ice  of  each  thickness. 

The  mathematical  equation  describing  changes  in  the  ice  thickness  distribution 
represents  mass  balance.  Ice  thickness  distribution  is  analogous  to  a 
set  of  equations  governing  changes  in  species  concentrations  in  a  mixture, 
except  that  the  ice  thickness  distribution  describes  an  infinite  set  of 
species  differentiated  only  by  their  thicknesses.  As  in  mixture  theory, 

species  interact  thermally  and  by  redistribution,  but  no  relative  motion 

is  allowed  between  species. 

The  constitutive  laws  are  introduced  to  relate  stress  in  the  ice 
pack  to  deformation.  This  stress  is  a  resultant  integrated  through 
thickness  much  like  the  stress  resultant  in  plate  or  shell  theory.  In 
addition,  the  formation  and  modification  of  leads  and  ridges  as  the  ice 
cover  deforms  changes  the  ice  thickness  distribution.  This  creation  of 
open  water  in  leads  and  the  formation  of  thick  ridges  from  thin  ice 
requires  a  constitutive  law  describing  the  production  and  conversion  of 
ice  due  to  mechanical  redistribution. 
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All  fundamental  aspects  of  this  mathematical  model  have  been  developed 
and  described  previously.  These  include  balance  of  mass,  momentum  and 
mechanical  energy.  Therefore,  we  present  the  governing  equations  without 
justifying  the  assumptions  made  to  obtain  each  balance  law.  The  constitutive 
laws,  however,  are  changed  enough  to  require  a  more  thorough  explanation. 

But  even  with  these  changes,  the  type  of  behavior  is  unchanged.  In  the 
remainder  of  this  section,  each  balance  law  and  constitutive  equation  is 
presented . 

3 . 2  MOMENTUM  BALANCE 

In  the  plane  of  motion  of  the  sea  ice,  the  momentum  balance  is 
expressed  as 

rav  =  -  mf  kXv  -  m^VH  +  7* a  ,  (3.1) 

where 

m  is  mass  per  unit  area  of  ice, 

v  is  ice  velocity  in  Che  horizontal  plane, 

• 

v  is  material  time  race  of  change  of  v  , 

Ta  is  traction  exerted  by  the  atmosphere  on  the  upper  surface  of  ice, 

T  is  traction  exerted  by  the  ocean  on  the  lower  surface  of  ice, 

~w 

f  is  the  Coriolis  parameter, 

g  is  acceleration  due  to  gravity, 

H  is  the  height  of  the  sea  surface, 
k  is  the  unit  vector  in  the  vertical  direction,  and 

O  is  the  Cauchy  stress  resultant  in  excess  of  hydrostatic  equilibrium 
(two-dimensional) . 

The  horizontal  position  vector  x  and  the  vertical  unit  vector 
k  are  expressed  in  a  right-handed  Cartesian  coordinate  system.  Velocity 
is  the  material  rate  of  change  of  position,  x  =  v  . 

The  acceleration  term  in  the  momentum  balance  is  generally  much 
smaller  than  the  other  terms  due  to  the  fact  that  the  acceleration  of 
the  ice  is  generally  small. 
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The  air  stress  is  a  dominant  driving  force  in  the  momentum  balance. 
For  a  typical  day,  it  is  the  largest  single  contribution  to  the  horizontal 
force  balance.  Because  air  stress  is  variable  from  day  to  day,  however, 
the  slowly  varying  ocean  current  can  often  provide  a  larger  contribution 
to  long-term  motions.  The  air  stress  x is  given  as  a  quadratic 
function  of  the  geostrophic  velocity  U  (Brown,  1976): 


T 

-a 


pc  |  U  i  B  U 
a  a  ~  ~a~ 


(3.2) 


where 


cos  a  -s in  & 


B  * 
-a 


sm  a  cos  a.. 


(3.3) 


and  p  *  1  kg/ra  is  air  density  and  c  =  0.0008  is  the  drag  coefficient. 
The  air  stress  is  turned  at  an  angle  a  =  28  counterclockwise  from  the 
applied  geostrophic  wind  U  (Albright,  1980;  Leavitt,  1980). 


The  effect  of  the  oceanic  boundary  layer  on  ice  motion  has  also 

been  modeled  as  a  quadratic  relation  between  the  traction  x^  exerted 

on  the  lower  surface  of  the  ice  and  the  ice  velocity  relative  to  the 

geostrophic  currents  v  (McPhee,  1975).  The  traction  is  applied  at  an 
o  ^ 

angle  8  =  23  to  the  relative  velocity: 


x 

-w 


pc  |v  -  v  I B(v  -  v  )  , 
w  w  -  ~g  1 - -g 


where 


(-cos  8 

-sin  8 


sin 


-cos 


(3.4) 


(3.5) 


3 

and  p  *  1000  kg/m  is  water  density  and  c  =  0.0055  is  the  drag  coefficient 
w  w 

Geostrophic  currents  are  considered  to  be  steady,  so  time  variations  do 
not  appear. 
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The  Coriolis  force  exerted  on  sea  ice  depends  on  the  ice  velocity, 
the  Coriolis  parameter  and  the  mass  of  sea  ice.  Within  the  AIDJEX 
model,  the  mass  of  the  sea  ice  is  determined  through  the  ice  thickness 
distribution  as  the  ice  density  (p  =  900  kg/m^)  times  the  mean  thickness. 

The  ocean  tilt  term  in  the  momentum  balance  arises  from  a  varying 
sea  surface  height  which, by  balancing  the  hydrostatic  head, determines 
the  geostrophic  flow  in  the  ocean: 

m^VH  =  -mf  kxv  .  (3.6) 

The  last  terra  of  the  momentum  balance  is  the  divergence  of  internal 
ice  stress.  The  stress  is  more  accurately  called  a  Cauchy  stress  resultant 
in  excess  of  hydrostatic  equilibrium  (Coon  et  al. ,  1974).  The  stress 
state  is  related  to  the  deformation  history  by  an  elastic-plastic  constitutive 
law  that  is  described  in  Section  3.5. 

3.3  MASS  BALANCE 

One  objective  of  AIDJEX  was  to  relate  the  observed  morphology  of 
the  ice  to  its  field  of  motion.  A  model  of  this  relation  has  been 
framed  in  the  theory  of  the  ice  thickness  distribution.  Here  the  ice 
pack  is  viewed  as  a  mixture  of  many  different  thicknesses.  A  balance 
equation  is  written  to  describe  how  each  thickness  category  changes  in 
response  to  the  deformation  and  the  thermodynamic  forcing.  While  this 
theory  does  not  attempt  to  describe  the  orientation  or  spacing  of  leads 
or  pressure  ridges,  it  does  account  for  their  net  effect  by  keeping 
track  of  the  amounts  of  open  water  and  of  thick,  pressure-ridged  ice. 

Many  of  the  large-scale  properties  of  the  ice  pack,  such  as  its  mass, 
its  strength  and  its  ability  to  insulate  the  atmosphere  from  the  ocean, 
depend  on  the  thickness  distribution.  It  can,  therefore,  be  viewed  as  a 
fundamental  state  variable. 

Balance  of  mass  is  the  principle  that  provides  an  equation  governing 
the  thickness  distribution.  Following  Coon  et  al.  (1974)  and  Thorndike 
et  al.  (1975)  we  write 
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(3.7) 


where 
G(x,  h,  t) 


fg(h,  o 


G  +  f  «  f  -  C  7*v 
g  3h 

is  the  fraction  of  area  in  the  vicinity  of  point  x  at 
the  time  t  that  is  covered  by  ice  thinner  than  the  ice 
thickness  h  , 

is  the  growth  rate  of  ice  and  may  be  different  for  each 
thickness  and  time,  and 

is  a  redistribution  function  that  describes  mechanical 
conversion  of  ice  between  categories. 


Changes  in  the  relative  area  covered  by  ice  in  each  category  are 

3G 

controlled  by  both  thermodynamic  growth  and  ablation  (f  -rr-)  and  by 

§►  oh 

mechanical  redistribution  (4')  .  The  remaining  term  (G  V*v)  is  needed 
to  account  for  changes  in  area  since  G  indicates  a  fraction  of  present, 
not  initial,  area.  The  time  rate  of  change  G  is  a  material  rate. 


The  fraction  of  open  water  in  an  area  is  described  by  the  zero¬ 
thickness  category  of  ice,  G(x,  0,  t)  ;  so  the  ice  compactness  is 
easily  found  by  subtracting  this  quantity  from  unity.  It  should  be 
noted  that  in  this  model  there  is  no  oceanic  heat  storage;  open  water 
formed  during  an  opening  or  shearing  episode  is  immediately  frozen  into 

thin  ice  during  the  "freezing  season"  when  f  (0,  t)  is  positive.  For 

S 

this  reason  and  because  it  is  difficult,  if  not  impossible,  to  differentiate 
between  open  water  and  thin  ice  using  available  remote  sensors,  we 
normally  do  not  try  to  separate  open  water  from  thin  ice,  say  less  than 
0.10  m  thick,  in  simulated  results  of  ice  behavior. 


The  rate  of  mechanical  production  of  ice  of  all  thicknesses  h  and 
thinner  is  ¥(h)  ,  where  we  suppress  dependence  on  x  and  t  which 
enter  implicitly  through  G  .  Ice  is  produced  in  this  range  of  categories 
by  the  formation  of  leads  and  by  conversion  from  other  categories  during 
ridging.  This  mechanical  redistribution  of  ice  between  categories 
changes  the  relative  fraction  of  area  covered  by  ice  of  each  thickness 
category  as  seen  in  Equation  (3.7).  The  amount  of  ridging  that  occurs 
at  a  location  depends  on  the  history  of  the  motion  of  that  ice.  The 
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thickness  distribution  serves  as  the  critical  measure  of  the  result  of 
that  history  of  deformation.  Therefore,  we  assume  that  the  instantaneous 
thickness  distribution  is  adequate  to  describe  the  history  and  that  only 
the  instantaneous  rate  of  deformation  has  an  influence  on  redistribution. 
Furthermore,  we  assume  that  processes  that  redistribute  ice  are  nonrecoverable 
and  thus  are  associated  with  permanent  deformations. 

3.4  MECHANICAL  ENERGY  BALANCE 

*  The  race  at  which  work  is  done  on  an  elemental  region  of  pack  ice 
by  all  forces  acting  on  the  ice  is  obtained  by  forming  an  inner  product 
of  the  ice  velocity  with  each  force  .n  the  momentum  balance  (Coon  and 
Pritchard,  1979;  Pritchard  et  al. ,  1979;  Thomas  and  Pritchard,  1980). 

As  a  result  of  this  mathematical  manipulation,  a  mechanical  energy 
balance  is  found  in  the  form 

;  +  Oe.Pi-po-pg  +  Pf-Pj  ,  U.6) 

where 

T 

<  =  1/2  mv  v  is  kinetic  energy  density  (the  superscript  T 

indicates  the  transpose), 

Ue  =  tr  ct  e  is  recoverable  strain  energy  density  (e  is  elastic 

strain) , 

T 

p  =  v  t  xs  power  input  to  the  ice  from  the  atmosphere  bv 

surface  traction, 

T 

p  “  ~v  t  is  power  loss  from  the  ice  to  the  ocean  bv  surface 

o  —  -»w 

traction  from  the  oceanic  mixed  layer, 

T 

Pg  «  mg  v  7h  is  the  rate  of  increase  of  gravitational  potential 

energy  due  to  ice  movement  up  the  sea  surface  slope, 

T 

Pf  *  v»(v  a)  is  divergence  of  the  stress  flux,  the  rate  at  which 

energy  is  transferred  horizontally  through  the  ice 
cover,  and 

P«  *  tr  0  D  is  dissipative  stress  power,  the  rate  at  which 

x.  -  ~p 

energy  is  dissipated  by  the  stress  state  a  during 

permanent  deformations  D 

~P 
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The  mechanical  energy  balance  relacionship  has  proven  to  be  an 
important  interpretive  tool.  It  serves  to  provide  a  measure  of  the 
level  of  activity  due  to  storms  and  to  understand  better  the  energy 
budget  of  polar  oceans.  Energetic  events  are  strongly  amplified  since 
is  proportional  to  the  wind  speed  cubed  when  ice  drift  is  a  linear 
function  of  wind  speed.  Furthermore,  the  small-scale  mechanisms  of 
energy  dissipation  occurring  during  ridging,  in  the  form  of  gravita¬ 
tional  potential  energy  increases  and  frictional  losses  from  sliding  of 
ice  blocks  into  ridges,  have  been  estimated  directly  (Rothrock,  1975). 

An  additional  shear  ridging  energy  sink  is  discussed  later  in  this 
chapter.  These  dissipative  mechanisms  are  related  to  the  large-scale 
dissipative  stress  power  to  constrain  the  constitutive  laws.  Thus,  if 
p  is  the  rate  of  dissipation  of  energy  from  all  subscale  mechanisms,  it 
must  equal  the  dissipation  rate  of  stress  during  deformation.  Therefore, 

PA  =  P  •  (3.9) 


3.5  CONSTITUTIVE  LAWS 

The  set  of  constitutive  laws  contains  relationships  describing 
stress  and  redistribution  in  terms  of  the  deformation  history.  Both  are 
intended  to  describe  the  large-scale  average  behavior  of  the  ice  cover 
over  length  scales  on  the  order  of  tens  of  kilometers.  This  average 
behavior  is  made  up  of  the  formation  of  many  randomly  oriented  cracks, 
leads  and  ridges,  assuming  enough  features  are  present  to  allow  a 
homogeneous,  isotropic  behavior.  The  material  behavior  is  assumed  to  be 
rate  independent,  as  determined  from  a  detailed  analysis  of  ridge 
formation  (Parmerter  and  Coon,  1972;  Coon,  1974).  The  material  response 
is  assumed  to  be  plastic  because  of  this  rate  independence  and  because 
of  the  fact  that  ridges  appear  to  have  limited  heights  when  formed  from 
an  ice  sheet  of  fixed  thickness.  The  material  model  to  be  described  has 
five  separate  elements  that  are  required  to  specify  completely  the 
internal  ice  stress  for  a  given  deformation  history.  These  elements, 
which  have  much  in  common  with  other  elastic-plastic  material  models, 
are:  yield  surface,  flow  rule,  elastic  response,  kinematic  relationship 
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and  hardening  law.  The  first  four  of  these  elements  are  now  discussed, 
while  the  fifth  is  introduced  after  describing  the  ice  redistribution 
law  on  which  it  depends. 


Figure  3.1.  Diamond  Yield  Surface  for  Large-Scale  Sea  Ice  Model 

The  stress  state  in  any  plastic  model  is  constrained  to  lie  within 
the  yield  surface.  For  an  isotropic  model,  this  function  depends  only 
on  the  stress  invariants  and  not  on  the  principal  direction.  Thus, 

♦  (Oj,  air  P*)  <  0  ,  (3. 10) 

u 

where  cT  =  1/2  tr  0  (negative  pressure),  a_  =  (1/2  tr  c'a')  (maximum 
l  -  *  ii  -  - 

shear),  0*0-  0^1  and  p  is  the  isotropic  compressive  strength. 

This  strength  will  be  shown  to  be  a  function  of  thickness  distribution. 

The  yield  surface  chosen  for  sea  ice  is  shown  in  Figure  3.1.  Along  the 
straight-line  portion  of  $  *  0  passing  through  the  origin,  the  stres- 
state  is  that  of  uniaxial  compression.  Simulations  of  sea  ice  dynamics 
during  January  1976  and  of  ice  flow  through  the  Bering  Strait  show  that 
it  is  important  to  use  this  "tensile  cutoff"  line  to  define  part  of  the 
yield  surface  (Pritchard,  1978;  F.eimer  et  al.  ,  1979;  Pritchard  and 
Reiraer,  1979).  The  other  straight  line  used  to  complete  the  yield 
surface  is  chosen  for  simplicity.  It  reduces  the  shear  strength  by  a 
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factor  of  2  compared  to  the  triangular  yield  surface  which  has  the 
largest  possible  shear  strength  compatible  with  tensile  cutoff,  a  concave 

: k 

surface  and  a  compressive  strength  of  p  ,  Since  the  effective  strength 
of -the  ice  cover  over  length  scales  of  tens  of  kilometers  cannot  be 

measured  by  direct,  controlled  experiments,  we  have  estimated  the  values 

* 

of  p  from  simulations  of  complicated  and  uncontrolled  large-scale  ice 

behavior.  A  careful  choice  of  wind,  current  and  ice  conditions  has 

allowed  simulations  to  be  performed  that  show  strengths  in  the  range  of 
4  5 

10  to  10  N/m  for  ice  conditions  common  to  the  Beaufort  Sea  during  the 

4 

winter.  It  is  worth  noting  that  strengths  below  10  N/m  have  little 
influence  on  ice  motions,  although  we  do  not  know  the  effect  of  ice 
stress  on  deformation  at  this  strength  level.  This  knowledge,  sparse  as 
it  is,  provides  one  test  of  the  reasonableness  of  the  strength  values 
determined  in  the  hardening/softening  plastic  model  to  be  presented  in 
the  remainder  of  this  section. 

The  next  three  elements  of  the  plastic  response  are  unchanged  from 
previous  forms  of  the  model  (Coon  et  al.,  1974;  Pritchard,  1975;  Coon, 

1980).  When  the  stress  state  in  the  ice  lies  inside  the  yield  surface, 

Chen  the  stress  a  is  a  function  of  the  elastic  strain  e  : 

a  =  (M1  -  M2)  1  tr  e  +  2M2e  ,  (3.11) 

where  and  M,  are  elastic  moduli  that  may  depend  on  the  thickness 

distribution.  The  elastic  moduli  are  normally  chosen  to  be  constant  and 
large  enough  so  that  elastic  response  is  a  representation  of  rigid 
behavior.  The  elastic  strain  satisfies  the  kinematic  relation 

e-We+eW=D-D  ,  (3.12) 

~P 

IT  IT 

where  D  »  (L  +  L  )  is  stretching,  W  *  —  (L  -  L  )  is  spin  and 

L  »  Vv  is  the  velocity  gradient. 

When  the  stress  state  is  on  the  yield  surface  <t>  =  0  ,  plastic 
stretching  occurs.  During  plastic  flow,  the  stress  is  constrained  to 
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Che  loading  surface.  The  plastic  screeching  D  satisfies  the  associated 

~P 

flow  rule 

5?  =  X  ,  where  X  >  0  (3.  13) 

which  requires  that  plastic  stretching  be  orthogonal  to  the  loading 
function  at  the  instantaneous  stress  state. 


The  redistribution  function  4*  retains  the  general  fora  developed 
by  Thorndike  et  al.  (1975).  This  constitutive  law  is  expressed  as  a 
function  of  permanent  or  plastic  stretching  D  in  the  form 


•a  (6)  +  a  (0)  W 
o  r  r 


(3.14) 


Invariants  of  the  irreversible  stretching  tensor  are  the  dilating 


Dj  =  tr  Dp  and  shearing  D^ 


(2  tr  D 'D ' )  ,  where  D*  =  D  -  1/2  DT1  . 

~p~p  ~p  ~p  1~ 


These  invariants  are  equal  to  the  sum  and  difference  of  the  principal 

stretching  components,  respectively.  Another  set  of  stretching  invariants. 

used  by  Nye  (1976),  is  more  valuable  for  describing  stretching  dependence 

?  2  U 

of  the  redistribution  function.  This  set  includes  D  »  (D“  +  DTT)~ 

p  III 

and  9  =  atan  (D^/Dj)  ,  which  are  measures  of  the  rate  of  stretching 

Dp  and  the  amount  of  shearing  relative  to  dilating  (0)  .  The  linear 

dependence  of  4/  on  Dp  ensures  that  the  redistribution  process  is 

rate  independent,  a  property  that  is  consistent  with  the  elastic-plastic 

stress  constitutive  law.  The  term  a  +  a  W  describes  the  production 

o  r  r 

of  ice  in  all  thickness  categories  thinner  than  h  per  unit  of  strain. 

The  coefficient  ct^  is  the  amount  of  open  water  produced  per  unit  of 
strain. 


Observations  of  the  redistribution  of  the  Beaufort  Sea  ice  cover 
indicate  that  a  reasonable  form  of  a  (6)  is  that  presented  in  Figure  3.2 
(Rothrock  and  Hall,  1975).  Therefore,  open  water  is  produced  not  only 
during  opening  deformations  (0  <_  9  <  tt/2)  but  also  during  shearing 
deformations  that  close  up  the  ice  slightly  (tt/2  <  0  <_  2ir/3)  .  The 
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Ridging,  af  Opening, 
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function  a  is  monoconic  decreasing,  indicating  that  more  open  water 
o 

is  formed  during  opening  than  during  closing  deformations,  a  result 
obviously  compatible  with  physical  intuition. 

The  production  of  ice  per  unit  strain  in  all  categories  thinner 
than  h  is  a^W  ,  which  describes  the  production  of  ice  from  ridging. 

The  coefficient  a^(9)  ,  like  aQ(6)  »  depends  on  the  relative  amounts 
of  shearing,  opening  and  closing.  Several  properties  of  are  important: 

a  is  positive  because  unridging  cannot  occur;  is  monotonic  increasing, 

indicating  that  more  ridging  occurs  during  closing  than  during  opening 
deformations;  and  finally, 

a  (9)  -  a  (9)  =  cos  6  (3.15) 

o  r 

to  satisfy  conservation  of  area. 

Following  Thorndike  et  al.  (1975),  ridging  is  described  by 

00 

n(h2)  =  J  Y(h1,  h2)  a(h^  dl^  ,  (3.16) 

0 

where  a(h^)  is  the  fraction  of  the  area  covered  by  ice  of  thickness 
h^  per  unit  thickness  that  participates  in  the  ridging  process  and 
n(h2)  is  the  fraction  of  the  area  of  ice  of  thickness  h2  per  unit 
thickness  that  is  produced  by  the  ridging  process.  The  function  Y(h,,  h,)  , 
called  the  ridging  process,  is  the  distribution  of  thicknesses  produced 
by  ridging  of  a  unit  area  of  ice  of  thickness  h^  .  The  infinite  limit 
of  integration  is  strictly  formal  in  that  the  integral  is  always  zero  at 
all  thicknesses  h^  >  h^^  ,  the  thickest  ice  present.  The  ridging  mode 
is  formed  by  combining  production  of  all  thinner  categories  while 
accounting  for  loss  of  ice  that  participates  in  the  ridging  process. 
Therefore, 

h 

Wr(h)  =  J  £n(h2)  -  a(h2)J  dh 
0 
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where  W  («°)  is  introduced  to  conserve  total  area.  Additional  important 
general  properties  of  the  redistribution  function  may  be  found  in  Thorndike 
et  al.  (1975).  Specific  choices  used  in  the  simulation  are  now  presented. 

The  ridging  process  is  chosen  to  reflect  the  assumption  that  all 
ridged  ice  is  k  times  its  original  thickness.  Thus, 

Y(hl*  h2}  =  k  S(h2  '  khl}  ‘  (3.18) 

Although  ice  produced  by  ridging  surely  is  in  a  distribution  of  categories, 
the  assumption  that  ice  of  thickness  h^  is  always  ridged  into  ice  of 
thickness  h^  =  kh^  is  attractive  for  its  simplicity.  The  value 
k  *  15  is  used  based  on  the  results  of  Parmerter  and  Coon  (1972).  The 
AIDJEX  model  (Coon  et  al. ,  1974)  used  k  =  5  ,  which  is  now  thought  to 
be  too  low.  Different  forms  of  ridging  processes  have  been  assumed  by 
Bugden  (1979)  and  more  recently  by  Hibler  (to  appear).  There  are  few 
data  available  to  test  these  assumptions. 

The  participation  function  a(h)  describes  the  fraction  of  area 
covered  by  ice  of  thickness  h  that  is  converted  by  the  ridging  process. 
This  function  is  bounded  above  by  the  fraction  of  area  covered  by  ice  of 
thickness  h  that  is  present.  Therefore,  it  is  useful  to  assume  that 

a(h)  =  b  rG(h)'  g(h)  ,  (3.19) 

-  J 

3G 

where  g  *  and  g(h)  dh  is  the  fraction  of  area  covered  by  ice 
thinner  than  h  +  dh  and  at  least  as  thick  as  h  (Thorndike  et  al. ,  1975). 
The  function  b  is  the  fraction  of  ice  present  that  participates  in  the 
ridging  process.  The  function  b  (Figure  3.3)  is  assumed  zero  for  the 
thinnest  GQ-th  percentile,  which  is  assumed  not  to  be  available  for 
ridging.  In  a  sense  we  protect  this  thinnest  ice,  assuming  that  it  lies 
in  the  openings  between  irregular,  large,  thick  floes  and  cannot  be 
ridged  without  fracturing  the  thicker  ice.  This  has  the  effect  of 
shifting  curve  b  to  the  right  in  Figure  3.3.  It  is  also  assumed  that 
the  function  b  decreases  with  increasing  values  of  G  .  That  is,  less 
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of  the  thicker  ice  is  available  for  ridging.  For  simplicity,  the  decrease 

is  linear  up  to  the  G^-th  percentile  and  then  constant  up  to  the  G^-th 

percentile.  The  small  amount  of  thicker  ice  available  for  ridging 

(G,  <  G  £  G^)  acknowledges  the  fact  that  some  thick  ice  is  fracturec 

and  piled  into  ridges.  In  Figure  3.3,  constants  have  been  chosen  as 

G  “  0.05,  G0  *  0.20,  G.  =  0.55,  and  B  =  0.90. 
o  Z  *  o 


Gi 

Cumulative  Thickness  Distribution,  G 


Figure  3.3.  Fraction  of  Ice  in  Each  Thickness  Category  Available  for  Ridging 


The  rate  at  which  energy  is  dissipated  by  the  deforming  ice  cover 

was  originally  related  to  the  changes  in  gravitational  potential  energy 

caused  by  mechanical  deformation  p  and  to  the  frictional  loss  in 

P 

ridging  p^  .  In  the  modified  formulation  we  introduce  another  energy 
sink  ps  to  account  for  energy  dissipated  by  shearing  deformations  that 
do  not  alter  the  ice  thickness  distribution. 

Analysis  of  each  of  the  first  two  energy  sinks  provides  a  similar 
relationship.  From  the  gravitational  potential  energy  changes  (Rothrock, 
1975), 


CD 


0 


(3.20) 
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where  c  »  4-  p g  and  p  =  (p  -  p)  .  By  substituting  the  assumed 
p  “  p  w 

w 

redistribution  4*  ,  this  becomes 


p  *  D  a  p 
P  P  r  P 


where 


p  =  c 
P  P 


/2  3W 

h  3irdh  • 


(3.21) 


(3.22) 


The  loss  due  to  friction  has  been  analyzed  and  estimated,  and  the  functional 
form  is  the  same  as  for  potential  energy  changes  (Rothrock,  1975)  : 


pf  =  DparPf  ’ 


where 


Pf  =  cf 


/  t4 


(h) 


(1/k) 


dh 


(3.23) 


(3.24) 


and 


9  r 

2  tan  p' 


P (k  -  1) 

Pw 


(3.25) 


where  u'  is  the  coefficient  of  sliding  friction  and  the  angle  p '  is 
the  slope  of  keels  of  pressure  ridges. 


To  determine  the  material  constants  c  and  c,  ,  mass  densities  and 

2  P  f 

the  gravitational  constant  g  =  9.80  m/s  are  used.  The  coefficient  of 
sliding  friction  is  u'  =0.35  ,  a  slight  increase  from  early  AIDJEX 
model  calculations.  The  increase  is  an  attempt  to  increase  yield  strength 
to  values  found  to  be  necessary  for  simulating  observed  ice  motions 
accurately.  This  value  is  within  the  range  of  reasonable  values,  but 
there  are  no  data  to  support  the  choice  of  any  particular  value  in  the 
ridging  process.  We  have  retained  the  value  of  about  39°  for  the  angle 
p’  ;  therefore,  tan  p'  =  0.8 


The  contributions  from  p  and  p^  are  combined  in  the  one  expression 


p  +  p,  =  D  a  p 
P  Kf  p  rK 


(3.26) 


where  p  =  p  +  p^  may  be  written  as 


*  * 

p  -  c 


/ 


h  a(h)  dh 


(3.27) 
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and 


c  =  4-  pgk  1  +  — -  -  , 
2  can  <J>  * 


£_ 


(k  -  1) 


For  the  chosen  material  constants,  we  find  that 


c*  =  43  x  103  N/m3 


(3.28) 


(3.29) 


The  energy  dissipated  by  friction  is  about  5.5  times  as  large  as  the 
energy  dissipated  by  potential  energy  changes.  Early  AIDJEX  model 
calculations  neglected  friction  and  had  a  very  weak  strength  as  a  result. 


The  final  energy  sink  has  not  been  analyzed  to  determine  how  shearing 
dissipates  energy.  We  assume  the  process  to  be  rate  independent  so  that 


p  is  linear  in 
s 


0^  ,  the  rate  of  plastic  deformation,  and  write 


Ps  -  Dpq  ,  (3.30) 

where  q  is  the  rate  of  dissipation  of  energy  per  unit  deformation 

rate.  It  may  depend  on  the  ratio  of  shearing  to  dilating,  the  thickness 

distribution  and  the  stress  state.  By  considering  the  relative  motion 

of  ice  on  either  side  of  a  velocity  discontinuity,  it  is  easy  to  imagine 

that  q  varies  as  the  product  of  the  ice  pressure  (-c^)  and  a  coefficient 

a  (0)  .  Then, 
s 


q  -  -Oj.  a,  ,  (3.31) 

but  we  do  not  make  a  direct  use  of  this  assumption.  Rather,  we  find  q 

as  a  residual  in  the  energy  balance.  The  requirement  that  q  be  nonnegative 

introduces  a  weak  restriction  on  our  choice  of  yield  surface  shape. 


When  we  equate  stress  power  to  rate  of  energy  dissipation  so  that 


P£  -  Pp  +  Pf  +  Ps  .  (3-32) 

this  mechanical  energy  balance  equation  provides  a  definition  of  the 
shearing  dissipation  q  for  all  plastic  stress  and  stretching  states. 

We  find  that 
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The  shearing  dissipation  depends  only  on  9  for  a  given  yield  surface 
and  flow  rule.  For  the  triangular  yield  surface  and  normal  flow  rule  as 
well  as  the  chosen  ridging  coefficient  ,  the  shearing  dissipation  is 
shqvm  in  Figure  3.4.  The  shearing  sink  is  never  negative,  a  condition 
necessary  to  prohibit  energy  generation.  It  is  seen  that  energy  is 

dissipated  in  shearing  when  tt/4  <  9  <  3ir/4  and  zero  outside  this 

* 

range.  The  value  of  q/p  has  a  maximum  of  0.25.  The  shearing  dissipation 
is  an  important  addition  to  the  constitutive  law.  With  the  addition  of 
a  shear  dissipation  sink,  less  open  water  is  formed  during  shearing,  and 
this  fact  prevents  or  reduces  softening  of  the  model  behavior  in  many 
deformations.  This  results  in  larger  values  of  strength  in  comparison 
with  values  produced  by  the  early  AIDJEX  model. 


Figure  3.4.  Relative  Amount  of  Shearing  Dissipation  Without  Redistribution 
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Chapter  4 

INPUT  DATA  REQUIREMENTS 


4 . 1  FREE-DRIFT 

At  each  location  at  which  ice  velocity  is  to  be  calculated,  the 
free-drift  model  requires  that  the  geostrophic  wind  velocity  and  deep 
ocean  current  velocity  be  specified.  The  geostrophic  wind  can  be 
replaced  by  the  surface  wind  if  the  drag  coefficient  c  and  turning 

3 

angle  a  are  changed  appropriately.  In  addition,  the  mass  area  density 
of  the  ice  cover  m  is  needed.  Daily  average  values  of  wind  and  current 
velocities  are  appropriate  so  that  ice  motion  trajectories  are  built  up 
by  vector  summation  of  day-to-day  displacements. 

An  individual  oil  spill  trajectory  may  be  calculated  by  inputting 
wind  and  current  vectors  at  the  location  of  the  spill  each  day.  However, 
if  a  general  capability  is  needed  for  following  any  oil  spill  in  a 
region,  e.g.,  the  southern  Beaufort  Sea,  then  the  wind  and  current 
fields  are  needed  everywhere  in  the  region.  To  date,  the  only  ocean 
current  information  available  is  the  steady  geostrophic  ocean  current 
field  determined  from  the  dynamic  sea  surface  topography. 

When  a  prediction  is  required,  winds  and  currents  must  also  be 
predicted  out  to  the  verifying  time.  The  best  geostrophic  wind  predictions 
are  useful  out  to  a  maximum  of  about  5  days.  When  model  simulations  are 
made  in  a  hindcast  mode  for  the  sake  of  determining  the  present  location 
of  an  oil  spill  that  cannot  be  observed  directly,  then  this  predictive 
time  limitation  does  not  apply. 

The  free-drift  velocity  is  relatively  insensitive  to  the  mass  area 
density  (or  mass  density  p  times  mean  thickness  h  ) .  Therefore,  it 
is  adequate  to  consider  only  climatological  changes  in  it.  In  this 
case,  no  data  are  required. 
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4.2  COMPLETE  ICE  DYNAMICS  MODEL 


Since  the  internal  ice  stress  transmits  forces  horizontally  through 
thq.  ice  cover,  ice  velocity  at  one  location  can  be  sensitive  to  driving 
forces  at  another  location.  Therefore,  the  complete  ice  dynamics  model 
requires  that  geostrophic  wind  and  current  fields  be  input  over  the 
entire  region  of  the  calculation.  These  fields  (the  geostrophic  ocean 
current  is  steady  state)  can  vary  continuously  in  time  to  allow  smooth 
transitions  during  numerical  integration,  but  time  histories  should  be 
filtered  to  remove  all  frequencies  with  periods  shorter  than  1  day. 
Barometric  pressure  fields  at  synoptic  intervals  (6  hours)  fill  this 
need  adequately. 

Calculating  ice  motions  within  a  region  requires  that  either  ice 
velocity,  traction  (the  traction  vector  describes  the  normal  and  shear 
force  applied  to  the  ice  edge)  or  a  combination  of  the  two  vectors  be 
specified  everywhere  around  the  boundary. 

At  those  boundaries  that  lie  along  shorelines,  the  ice  motion  is 
assumed  to  be  zero.  Slippage  of  the  ice  along  the  shoreline  is  auto¬ 
matically  calculated  by  the  model  as  a  plastic  flow  within  the  ice 
nearest  shore.  Along  those  boundaries  corresponding  to  a  free  ice  edge, 
the  traction  can  be  prescribed  to  be  zero.  Since  the  modex  calculates 
motion  of  the  ice,  this  ice  edge  is  automatically  followed.  Neither  of 
the  above  boundary  conditions  requires  input  of  observations. 

When  the  region  of  interest  is  bounded  by  a  line  lying  within  the 
interior  of  the  ice  cover,  either  the  velocity  or  traction  along  the 
line  is  required.  For  a  hindcast  calculation,  data  buoys  can  be  deployed 
to  provide  such  input  data.  This  method  has  been  used  in  the  past  (e.g., 
Pritchard,  1981).  Another  possible  aproach  is  to  extend  the  region  of 
interest  to  a  coastal  boundary  or  a  free  ice  edge.  This  solution, 
however,  increases  the  requirement  for  additional  wind  and  current  data 
and  can  introduce  substantial  errors  if  these  driving  forces  are  known 
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only  poorly.  Although  Pritchard  and  Thomas  (1980)  show  that  errors  in 
calculated  ice  velocity  decay  rather  strongly  with  distance  from  the 
boundary,  the  question  of  how  to  specify  conditions  at  a  boundary  within 
the  ice  cover  remains  an  open  research  question  at  this  time. 

For  a  forecast  calculation,  the  boundary  velocity  must  be  predicted. 
One  approach  in  predicting  the  boundary  velocity  is  to  assume  persistence 
and  to  input  the  last  known  values  after  forward  extrapolation.  Free- 
drift  predictions  could  also  be  used  as  boundary  velocities,  but  the 
level  of  performance  of  the  complete  model  then  is  degraded  and  at  times 
is  nc?  better  than  the  free-drift  model  (Coon  et  al.,  1977). 

Finally,  initial  conditions  require  that  the  thickness  distribution, 

velocity  and  stress  be  specified  everywhere.  The  latter  two  variables 

have  little  long- lasting  effect  on  response  and  can  be  input  as  nominal 

values,  say  v  =0,0  =  -  m  1  ,  where  v  and  a  are  initial  values 

~o  -  ~o  l  ~  ~o  ~o 

of  v  and  a  .  But  the  initial  thickness  distribution  affects  behavior 
for  several  weeks  typically.  The  thicker  ice  concentrations,  say  thicker 
than  1.5  m,  can  be  input  from  climatological  values.  The  thin  ice 
requires  that  satellite  or  aircraft  imagery  be  obtained  to  determine  the 
ice  compactness  and  the  concentrations  of  thin  ice  in  several  categories. 
A  combination  of  visual,  infrared,  radar  and  microwave  sensors  is  useful 
to  provide  all  weather  coverage  and  to  identify  the  various  categories 
needed. 
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Chapter  5 
MODEL  PERFORMANCE 


5 . L  DEFINITION  OF  ERRORS 

An  evaluation  of  the  free-drift  and  complete  ice  dynamics  models 
requires  some  test  of  their  accuracy  in  predicting  observed  motions 
(Thorndike  and  Cheung,  1977).  Output  from  these  models  is  in  the  form 
of  daily  velocity  vectors.  Velocities  are  integrated  to  produce  trajectories. 
For  comparison,  the  actual  observed  trajectories  are  differentiated  to 
produce  daily  velocity  vectors.  The  vector  difference  between  these  and 
the  model  velocities  is  the  error  measure  used. 

Daily  error  vectors  exhibit  a  good  deal  of  variability.  In  order 
to  predict  the  error  expected  in  using  one  of  the  models,  variance 
ellipses  have  been  constructed.  These  give  an  area  within  which  the 
error  vector  is  expected  to  lie,  at  some  confidence  level,  for  each 
geographic  location  where  a  time-series  comparison  is  available. 


Computation  of  the  mean  error  and  standard  deviation  of  the  error 
follows  standard  statistical  definitions.  We  let  v  be  the  calculated 
daily  velocity  and  w  be  the  observed  daily  velocity.  Although  no  sub- 
or  superscript  is  used,  each  variable  is  evaluated  daily.  The  daily 
velocity  error  is 

£  =  v  -  w  (5.1) 

During  1  month,  the  mean  error  is 

S-s2  5  '  (5-2> 


where  N  is  the  number  of  daily  observations  during  the  month.  The 
standard  deviation  of  the  error  is 


(5.3) 
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Because  of  large  errors  in  current  velocities  v  ,  it  is  useful  to 

"**8 

isolate  the  effect  on  total  daily  error  z  •  This  separation  is  useful 

when  different  error  measures  represent  errors  from  different  physical 

processes.  To  this  end,  we  introduce  the  relative  ice  velocity  G  by 

separating  the  geostrophic  current  v  : 

g 


G  +  v 
-g 


(5.4) 


Then  the  error,  mean  error  and  standard  deviation  actually  computed  are 


e  =  G  -  w 


-  1  V 

e  =  —  '  e 
~  N  Z-  ~ 


(5.5) 

(5.6) 

(5.7) 


Combining  Equations  (5.1),  (5.4)  and  (5.5),  we  obtain 


£  =  e  +  v 
~g 


(5.8) 


and  it  follows  that  when  geostrophic  ocean  currents  are  included  in 
calculating  the  ice  drift,  the  mean  daily  error  is 


£  =  e  +  v 


(5.9) 


and  the  standard  deviation  of  the  error  e  +  v  is 

~g 


-( 


+2  cov  ( 

e  v 

g 


!•  V) 


(5.10) 


where  the  covariance  is  bounded  as 


-S  S  <  cov  (e,  v  )  <  S  S 
e  v  —  ~  ~g  —  e  v 

g  g 


(5.11) 
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When  ocean  currents  are  constant  along  the  month-long  trajectory,  or 
nearly  constant,  then  satisfies 


v  =  v  and  S  =0 
'S  'S  "g 


(5.12) 


The  standard  deviation  of  daily  errors  £  becomes 


S  =  S 
£  e 


(5.13) 


If  does  vary,  the  standard  deviation  of  the  daily  error,  using 

Equations  (5.10)  and  (5.11),  is  found  to  be  bounded  by 


s  <  s  -  s 

£  —  e  v 

g 


(5.14) 


which  shows  the  range  of  influence  that  ocean  current  variations  may 
have  on  the  standard  deviation  of  errors  S 

£ 

Variance  ellipses  are  constructed  as  follows.  Rectangular  coordinates 
are  used  to  avoid  difficulties  inherent  in  working  with  azimuths.  The 
variance-covariance  matrix  at  each  point  for  each  month  is  also  computed: 


S2  E 
x  xy 


S  S 
xy  y 


where 


2  IV  —  2 

s,  •  fi:  <*-«> 


s5  •  srrl  ^ 


(5.15) 


sxv  ■  3TT  I  <--»>  O-y) 


A  bivariate  normal  distribution  has  the  probability  density  function 


f(x,y) 


(5.16) 


2 irS  S  /  1-r 
x  y  V 


where 


Q  = 


(1-r2) 


(x-x)2 


2r (x-x) (y-y)  (y-y)2 

S  S  _2 

x  y  sy 


(5.1" 


wich  the  mean  (x,  y)  and  correlation  r  =  S 


xy 


The  locus  of  Q  =  c  (c  is  constant)  defines  an  equiprobability 

ellipse.  The  probability  p  that  a  random  point  on  the  x-y  plane 
will  fall  within  the  ellipse  is 


P  = 


(5.18) 


The  probability  that  the  magnitude  of  the  error  vector  lies  within  the 
range  defined  by  the  standard  deviation  is  p  *  0.683.  This  corresponds 
to  a  value  of  c  =  1.515.  Values  of  c  for  p  =  0.50,  0.90  and  0.99 
are  respectively  1.177,  2.146  and  3.035.  The  variance  ellipses  plotted 
represent  one  standard  deviation  or  p  =  0.683.  Values  of  c  determine 
the  size  of  the  ellipse,  so  the  given  values  may  be  used  to  scale  ellipses 
at  different  probability  levels. 

A  transformation  of  variables  is  made  to  simplify  locating  the 
equiprobability  ellipses.  Morrison  (1976),  among  others,  has  described 
the  procedure.  First,  a  translation  is  made  and  a  new  coordinate  system 
is  defined  as 


x  =  x  -  x  ,  (5 . 19) 

~  ~o  ~o 

where  xq  refers  to  the  original  coordinates  and  x  now  refers  to  the 
new  coordinates.  The  mean  in  the  new  coordinate  system  is  now 
(x,  y)  *  (0,  0)  . 

The  variance-covariance  matrix  S  is  not  affected  by  a  linear 
transformation.  The  determinant  equation 

|S  -Al|  =  0  (5.20) 
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is  solved  for  Che  characteristic  roots  (or  eigenvalues)  A^  and  A^ 
and  the  equation 

(S  -  AI)e  =  0 


(5.21) 


is  solved  for  the  characteristic  vectors  (or  eigenvectors)  e‘  =  (cos  a,  sin  a) 

T  ~  * 

and  =  [cos  (a-Ki/2) ,  sin  (a+ir/2)]  where  a  is  the  angle  the  first 

principal  axis  makes  with  the  x-axis.  A  rotation  of  coordinates  gives 


T 

u  =  x  e , 


v  =  x  e^ 


(5.22) 


where  the  new  vector  components  u  and  v  are  independent  and  un¬ 
correlated  with  the  variances 


2  , 
a  =  A. 
u  1 


(5.23) 


a2  =  A 

v  2 


In  this  last  coordinate  system.  Equation  (5. 17)  becomes 

2  2 

Q  -  y~  +  y- 
AI  A2 

and  when  the  substitution  Q  =  c^  is  made  we  have 


(5.24) 


2.  2. 

c  A  ^  c  A  2 


=  1  , 


(5.25) 


which  is  the  equation  of  an  ellipse  with  major  and  minor  axes  of  lengths 
2c  J A^  and  2cJ  A2  ,  respectively. 
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The  entire  statistical  procedure  is  coded  in  Fortran.  Results 
from  the  various  modeling  runs  are  read  from  tapes  and  compared  with 
observed  trajectories  on  data  cards.  Positions  of  data  stations  do  not 
always  correspond  to  grid  nodes  in  the  model.  For  these  cases,  the  node 
closest  to  the  station  is  used  as  if  they  are  the  sane  point.  Output 
includes  trajectories  for  the  complete  ice  dynamics  model,  the  free- 
drift  model  and  observed  motion.  Average  velocity  error  vectors  are 
computed  from  the  difference  in  trajectories  at  the  end  of  the  period. 
The  semimajor  and  semiminor  axes  of  the  variance  ellipses  are  given  with 
their  orientation.  Standard  deviations  are  also  computed. 

5.2  RESULTS  OF  SIMULATIONS 

The  free-drift  and  complete  ice  dynamics  models  have  been  tested 
against  observed  ice  motions  during  four  periods.  During  the  AIDJEX 
experiment,  manned  camps  and  data  buoys  were  used  to  track  motion 
of  pack  ice  in  the  Beaufort  Sea.  Observations  were  made  from  April  1975 
to  May  1976  on  at  least  a  daily  basis.  From  these  observations,  four 
time  periods  are  chosen  to  test  the  models:  13-28  April  1975  (75RUN2B) , 
17-25  May  1975  (75RUN1F),  27  January  through  12  February  1976  (76RUN5C) 
and  June  through  December  1975.  A  free-drift  calculation  using  the  same 
input  data  is  made  for  each  run.  No  changes  are  made  in  the  free-drift 
calculation  from  run  to  run. 

5.2.1  75RUN2B:  13-28  April  1975 

Only  the  position  of  the  Big  Bear  camp  (BB)  is  available  for 

this  calculation.  The  modeled  area.  Figure  5.1,  extends  from  Banks 

Island  to  Point  Barrow  in  a  roughly  rectangular  area.  The  boundary  is 

constrained  to  be  motionless  alongshore  and  to  follow  free-drift  along 

the  outer  grid  nodes.  Ice  strength  is  determined  by  ice  thickness 

distribution.  An  initial  distribution  was  estimated  from  NCAA  satellite 

images.  As  the  ice  deforms,  the  thickness  distribution  changes  causing 

hardening  or  softening.  Thermal  growth  of  ice  during  the  period  of  the 

experiment  was  included  in  the  thickness  distribution.  The  ice  strength 
*  3  A 

varies  from  p  =  5  x  10  N/m  to  1  x  10  N/m  .  These  strengths 
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Figure  5.1.  Region  Modeled  in  75RUN2B.  Location  of  Big  Bear  (BB)  is 

Shown  on  11  April  1975.  [Shaded  Regions  Identified  by  Roman 
Numerals  Were  Used  by  Coon  et  al.  (1977)  to  Describe  Different 
Ice  Conditions.] 
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are  too  low  to  represent  large-scale  strength  accurately.  During  this 

time,  a  low  pressure  system  moves  eastward  across  the  Beaufort  Sea  at 

~200  km/day.  This  system  gives  rise  to  a  moderately  strong  air  stress 

2 

of.~0.2  Pa  (1  Pa  =  1  N/m  )  near  Big  Bear.  Further  details  of  this 
experiment  are  given  by  Coon  et  al.  (1977). 

Figure  5.2  shows  observed  and  simulated  trajectories  for  the  Big 
Bear  camp  from  13  through  28  April  1975.  The  actual  trajectory  is  shown 
along  with  a  trajectory  produced  by  the  free-drift  model  and  one  produced 
by  the  complete  ice  dynamics  model  (the  early  AIDJEX  model  is  used  for  the 
complete  model  in  this  simulation).  For  this  computer  run,  75RUN2B,  the 
strength  input  to  the  AIDJEX  model  is  so  low  that  it  behaves  very  similarly 
to  the  free-drift  model.  Both  models  track  the  actual  motion  very  well 
during  the  first  7  days  of  the  trajectory.  On  the  eighth  day,  the  camp 
slows  and  begins  an  erratic  motion  while  the  models  predict  continued 
motion  to  the  southeast.  During  the  first  3  days,  the  motion  is  offshore 
to  the  northeast.  Then  a  steady  southeast  motion  begins  which  tends  to 
pile  the  ice  into  the  shore  of  Banks  Island.  By  the  seventh  day  the 
effect  of  compressing  the  ice  against  shore  reaches  Big  Bear.  Tue 
strength  of  the  ice  is  sufficient  to  prevent  any  further  motion,  and 
the  wind  remains  more  or  less  steady.  The  free-drift  model  does  not 
account  for  ice  strength  and  thus  predicts  continued  motion.  The  strength 
in  this  version  of  the  complete  ice  dynamics  model  is  so  low  as  to  be 
insignificant  and  it  too  allows  continued  motion  to  the  southeast. 

These  trajectories  were  used  to  calculate  a  mean  velocity  error 
vector.  This  is  simply  the  error  in  the  trajectory  after  15  days  divided 
by  the  total  time.  Variance  ellipses  are  calculated,  see  Figure  5.2, 
for  the  variance  in  daily  velocity  error  vectors  about  the  mean  error 
vectors.  The  ellipses  represent  the  region  in  which  68.3  percent  of 
daily  velocity  error  vectors  are  expected  to  fall.  The  elongation  of 
the  variance  ellipses  is  probably  caused  by  the  prevalent  motion  in  the 
direction  of  the  elongation.  The  variance  is  equal  to  the  sum  of  the 
squares  of  the  major  and  minor  axes  of  the  ellipse.  The  standard  deviation 
is  given  along  with  daily  velocity  error  vectors  in  Table  5.1. 
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N  OBSERVED 


v<m/s) 


u(m/s) 


=d 


Free-OrrfT  Err  Of  Vector 
and  Variance  Ellipse 


Model  Enof  Vector 
and  Variance  Ellipse 


Figure  5.2.  Ice  Trajectory  at  Big  8ear  During  13-28  April  1975  (75RUN2B) 
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TABLE  5.1.  ERROR  STATISTICS  FOR  COMPLETE  ICE  DYNAMICS  AND  FREE-DRIFT  MODELS  (x,  y) 


*CAB  =>  Central  Arctic  Basin;  TRANS  =  Transition  Zone;  NS  =  Nearshore  Region. 


5.2.2  75RUN1F:  15-25  May  1975 

During  this  time,  positions  are  available  for  four  manned  camps  and 
six  data  buoys  surrounding  them  (Figure  5.3).  The  four  manned  camps  are 
designated  Caribou  (Ca) ,  Big  Bear  (BB) ,  Blue  Fox  (BF)  and  Snow  Bird  (SB). 

All  the  camps  are  on  the  ice  pack  of  the  central  Arctic  basin  over  500  km 
offshore.  Motions  of  the  data  buoys  are  used  to  provide  boundary  conditions 
for  the  AIDJEX  model  calculation  and  are  not  used  in  claculating  trajectory 
errors.  The  time  period  for  this  comparison  is  chosen  partly  because  of 
the  availability  of  Landsat  imagery,  which  provides  an  independent  check 
of  ice  motions.  A  comparison  of  Landsat  observed  and  modeled  motions  is 
given  by  Hall  (1980).  The  ice  is  modeled  as  a  perfectly  plastic  material 
with  a  constant  (no  hardening)  uniaxial  compressive  strength  of 

*  4 

p  =  4  x  10  N/m  .  Further  details  on  this  model  and  the  effect  of 
strength  are  given  by  Pritchard  (1978). 

There  are  four  manned  camp  trajectories  available  for  comparison 
with  75RUN1F  and  free-drift  model  calculations.  The  trajectories  are 
shown  for  the  four  manned  camps  in  Figures  5.4  through  5.7.  These  cover 
Che  period  from  17  to  25  May  1975.  The  motions  are  generally  to  the 
northwest  for  the  first  3  days  then  they  turn  slowly  to  the  southeast. 

This  southeast  drift  motion  of  the  camps  is  similar  to  that  for  the 
period  of  75RUN2B.  There  is,  however,  no  time  when  the  camps  stop 
moving  for  any  length  of  time.  This  is  reflected  in  the  error  and 
variance  plots. 

Both  models  track  the  camp  motion  fairly  well,  with  the  exception 
of  the  free-drift  trajectory  calculated  for  the  Blue  Fox  camp  (Figure  5.6). 
It  is  interesting  to  note  that  the  orientation  of  the  variance  ellipses 
does  not  vary  much  from  the  northwest-southeast  elongation.  The  magnitudes 
of  error  vectors  and  standard  deviations  are  given  in  Table  5.1.  Both 
the  mean  errors  and  standard  deviations  are  smaller  than  for  75RUN2B, 
indicating  that  strength  has  some  importance.  This  has  been  verified  by 
Hall  (1980)  when  comparing  results  of  other  simulation:  using  lower 
values  of  strength  (75RUN1E).  In  Hall’s  study,  the  constant  strength 
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Figure  5.3.  Region  Modeled  in  75RUN1F.  Locations  of  the  Caribou  (Ca),  Big  Bear  (BB). 

Blue  Fox  (BF)  and  Snow  Bird  (SB)  Camps  and  the  Six  Data  Buoys  (Circles  on 
Boundary)  Are  Shown  on  15  May  1975. 
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Figure  5.4.  Ice  Trajectory  at  Caribou  During  17-25  May  1975  (75RUN1F) 
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Figure  5.5.  Ice  Trajectory  at  Big  Bear  During  17-25  May  1975  (75RUN1F) 
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Figure  5.6.  Ice  Trajectory  at  Blue  Fox  During  17*25  May  1975(75RUN1F) 
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Figure  5.7.  Ice  Trajectory  at  Snow  Bird  During  17-25  May  1975  (75RUN1 F) 
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AIDJEX  model  calculation  is  compared  with  Landsat  data.  He  calculates  a 
"figure  of  merit"  from  velocities  for  1  day  at  10  positions  on  the  ice: 


(5.26) 


where  M  and  L  are  model  and  Landsat  velocity  vectors,  respectively. 
The  standard  deviation  is 


so 


(5.27) 


(5.28) 


and  for  10  cases,  N  =  10  ,  S£  =  1.054  A^  .  Hall  gives  A^  =  1.4  km/day 
corresponding  to  S^.  =  0.017  m/s  .  This  value  compares  well  with  our 
estimated  S^,  of  0.018  m/s  for  the  75RUN1F  complete  ice  dynamics  model 
calculation. 


5.2.3  76RUN5C:  27  January  Through  12  February  1976 


In  this  last  run,  the  complete  ice  dynamics  model  as  described  in 

Chapter  3  is  used.  Ice  motions  within  the  modeled  region  are  monitored 

by  three  manned  camps  and  eleven  data  buoys  scattered  throughout  the 

region  (see  Figure  5.8).  Boundary  conditions  are  fixed  along  the  shore 

from  Banks  Island  to  Point  Barrow  and  constrained  to  follow  the  motions 

of  five  data  buoys  located  around  the  outer  nodes  of  the  region.  As  in 

75RUN2B,  an  ice  thickness  distribution  is  used  to  determine  ice  strength. 

The  important  difference  is  that  the  initial  strength  is  two  orders  of 

*  5 

magnitude  higher,  p  -  10  N/m.  A  more  detailed  discussion  of  the  results 
that  is  provided  here  is  presented  by  Pritchard  (1981). 


The  period  of  this  simulation  is  chosen  because  large  temporal  and 
spatial  variations  of  air  stress  due  to  the  passage  of  two  storms  were 
thought  to  provide  a  good  test  of  the  model  under  extreme  conditions. 
Both  strong  onshore  (compressive)  and  offshore  (tensile)  motions  were 
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Figure  5.8. 


Region  Modeled  in  76RUN5C.  Locations  of  the  Caribou  (Ca), 
Snow  Bird  (SB)  and  Blue  Fox  (BF)  Camps  and  the  Data  Buoys 
(Small  Numbered  Circles)  Are  Shown  on  27  January  1976. 


observed  as  well  as  sharp  discontinuities  in  ice  motion  giving  rise  to 
shearing  in  the  stamukhi  zone.  High  winds  accompanying  storms  can  cause 
large  ice  motions,  because  of  the  quadratic  drag  law,  or  no  motion  at 
all,  when  the  ice  hardens. 

Fourteen  station  trajectories  scattered  throughout  the  Beaufort  Sea 
are  available  for  comparison  with  simulated  motions.  Rather  than  present 
trajectories  for  all  these  stations  here,  four  are  chosen  to  represent 
distinct  regions  in  the  pack  ice.  (Table  5.1,  however,  includes  the 
error  statistics  for  all  fourteen  trajectories.)  The  Caribou  manned 
camp  (Ca)  is  located  in  the  middle  of  the  Beaufort  Sea  well  into  the 
central  Arctic  basin  (see  Figure  5.8).  In  the  winter,  this  is  a  region 
of  well-consolidated,  uniform  pack  ice.  Station  17  is  a  RAMS  buoy 
located  off  the  west  coast  of  Banks  Island.  This  is  a  region  where 
large  leads  are  continually  opening  and  closing  in  response  to  the 
prevailing  winds,  which  alternately  drive  the  ice  pack  east  and  west. 
Station  20  is  a  RAMS  buoy  in  the  nearshore  region.  The  ice  in  this 
area  is  relatively  immobile  compared  to  the  rest  of  the  pack  because  of 
the  influence  of  the  shoreline.  The  last  region  considered  is  the 
highly  sheared  stamuKhi  zone.  This  is  the  transition  zone  between  the 
mobile  ice  pack  and  the  fast  ice  nearshore.  Station  16  is  used  to 
represent  the  motion  in  this  region. 

In  each  of  the  four  comparisons,  the  motion  predicted  by  free-drift 
is  much  larger  than  the  observed  motion.  This  is  seen  in  both  the 
trajectories  and  the  variance  ellipses  in  Figures  5.9  through  5.12.  The 
period  from  27  January  to  12  February  1976  includes  two  large  storms. 

The  first  of  these  drives  the  ice  to  the  west;  the  second  back  to  the 
east.  For  both  storms,  the  free-drift  calculations  result  in  motions 
chat  are  clearly  too  large.  The  resulting  variance  ellipses  are  quite 
large,  with  orientations  roughly  NNE  to  SSW  and  a  maximum  variance  of 
about  0.10  m/s.  The  complete  ice  dynamics  model  does  a  much  better  job 
of  predicting  motions  at  this  time  of  year.  This  is  not  surprising 
since  the  ice  has  considerable  strength  in  the  winter. 
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Figure  5.10.  Ice  Trajectory  at  Station  17  During  27  January  - 
12  February  1976  (76RUN5C) 
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Figure  5.11.  Ice  Trajectory  at  Station  20  During  27  January  • 
12  February  1976  (76RUN5C) 
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There  are  some  interesting  effects  of  geographic  location  on  the 
complete  ice  dynamics  model  errors.  The  standard  deviations  given  in 
Table  5.1  show  that  variances  are  consistently  smallest  in  the  central 
Arctic  basin,  largest  in  the  transition  zone  and  have  intermediate 
values  in  the  nearshore  region.  This  is  a  consequence  of  the  fixed 
boundary  conditions  at  the  shore.  The  model  does  best  furthest  removed 
from  boundary  constraints.  As  the  fixed  boundary  begins  to  influence 
motions,  the  model  does  poorly  because  of  Che  discontinuous  nature  of 
the  transition  region.  In  the  nearshore  area,  the  model  still  predicts 
some  motion  while  the  ice  in  fact  moves  hardly  at  all.  While  relative 
errors  in  this  area  may  be  large,  the  motions  are  small  and  the  error 
vectors  are  thus  not  as  large  as  in  the  transition  region.  The  transition 
region  between  the  fast  nearshore  ice  and  the  pack  ice  is  the  most 
difficult  in  which  to  predict  motions.  This  is  a  region  of  strong 
discontinuities  in  ice  motion.  In  the  stamukhi  zone,  shear  ridging  and 
shear  leads  are  confined  to  narrow  strips  which  are  not  resolved  by  the 
model.  The  region  west  of  Banks  Island  frequently  has  large  leads 
opening  and  closing,  giving  rise  to  large  discontinuities  in  ice  motion. 

It  is  interesting  to  note  that  the  variance  ellipses  for  these  regions 
are  consistently  eloi  zated  parallel  to  the  shoreline:  east-west  off 
Point  Barrow  and  north-south  off  Banks  Island  (except  for  Station  17 
which  could  be  due  to  its  location  with  respect  to  Amundsen  Gulf) . 

This  is  largely  because  large  onshore  motions  are  not  possible. 

5.2.4  June  Through  December  1975 


This  long-term  free-drift  calculation  covers  the  period  from  late 
spring,  when  the  ice  is  still  fairly  strong,  through  summer,  when  the 
ice  has  little  strength  (due  to  large  fractions  of  open  water)  and  free- 
drift  is  expected  to  work  well,  and  into  the  fall  when  the  ice  regains 
its  strength. 

The  observed  positions  of  three  manned  camps  are  averaged  and  used 
to  compute  monthly  ice  velocity  vectors.  These  are  compared  with  monthly 
velocities  from  the  free-drift  model  to  compute  error  magnitudes  and 
standard  deviations  (Thomas  and  Pritchard,  1979).  The  results  are  given 
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in  Table  5.2.  Error  magnitudes  are  fairly  uniform  from  June  through 
December  at  about  0.01  m/s.  The  standard  deviation  on  the  other  hand  is 
significantly  higher  in  November  and  December  (0.06  to  0.07  m/s)  than  in 
the  summer  (0.03  m/s). 

Table  5.2.  LARGE-SCALE,  MONTHLY,  FREE-DRIFT  ERROR  MAGNITUDES, 

JUNE  THROUGH  DECEMBER  1975.  (Thomas  and  Pritchard,  1979) 


Error 

June 

July 

Aug. 

Sept. 

.Oct. 

Nov. 

Dec. 

| e |  m/s 

0.006 

0.007 

0.009 

0.011 

0.006 

0.006 

0.006 

S  m/s 

0.025 

0.033 

0.034 

0.034 

0.043 

0.072 

0.066 

5.2.5  Twenty-Five  Year  Variability  of  Motions 

As  part  of  the  research  conducted  for  OCSEAP,  a  set  of  free-drift 
ice  motions  has  been  calculated  on  a  monthly  basis  (Thomas  and  Pritchard, 
1979) .  The  input  for  these  calculations  was  geostrophic  winds  based  on 
25  years  of  daily  sea  level  pressures  (1953-78)  provided  by  the  National 
Center  for  Atmospheric  Res jarch  (Jenne,  1975).  Free-drift  motions  were 
calculated  monthly  in  the  Beaufort  and  Chukchi  seas  by  accumulating  daily 
motions.  For  every  month,  these  data  were  then  averaged  over  the  25 
years  to  give  mean  monthly  motions  and  variance  ellipses. 

A  typical  variance  ellipse  has  major  and  minor  axes  of  0.06  m/s  and 
0.02  m/s  for  a  standard  deviation  of  about  0.06  m/s.  The  smallest 
values  occur  in  June  and  July  with  standard  deviations  as  low  as  0.03  m/s 
near  Prudhce  Bay.  The  largest  values  are  in  the  Chukchi  Sea  in  November- 
February  where  one  standard  deviation  may  reach  0.09  m/s. 
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In  the  Beaufort  Sea  in  May,  the  standard  deviation  of  free-drift 
ice  motions  was  '■0.04  m/s.  This  may  be  compared  with  a  standard  deviation 
of  errors  of  0.02  to  0.05  m/s  from  the  free-drift  calculations  during 
13-28  April  1975,  17-25  May  1975  and  the  June  through  October  portion 
of  the  June  through  December  1975  free-drift  calculations.  An  ice 
trajectory  for  1  month  calculated  from  the  25-year  averaged  ice  motions 
will  thus  be  about  as  accurate  as  a  complete  free-drift  calculation  for 
this  time. 

The  situation  is  similar  during  the  period  from  27  January  to 
12  February  1976.  The  standard  deviation  from  the  free-drift  calculation 
is  about  0.10  m/s  compared  to  expected  variations  from  year  to  year  of 
only  0.09  m/s. 

Only  the  complete  ice  dynamics  model  gives  a  significant  improvement 
over  the  25-year  average  motions;  the  standard  deviation  drops  to  0.036  m/s 
which  is  less  than  half  of  the  variability  in  ice  motions. 

5.3  DISCUSSION 

The  ice  drift  models  may  now  be  evaluated  in  terms  of  their  ability 
to  predict  daily  ice  trajectories.  The  accuracy  of  prediction  is 
necessarily  tied  to  the  accuracy  of  the  input  winds  and  currents.  Other 
errors  arise  when  the  model  is  inadequate,  as  when  free-drift  is  assumed 
in  winter.  With  the  statistics  assembled,  it  is  possible  to  estimate 
how  well  the  two  models  will  perform  at  different  times  of  year  and  what 
the  influence  of  wind  and  current  errors  will  be. 

5.3.1  Effect  of  Errors  in  Input  Winds  and  Currents 

Ultimately,  Coast  Guard  operations  will  require  meaningful  accumulated 
trajectory  errors  based  on  available  wind  and  current  data.  It  is  not 
possible  to  assess  errors  in  the  winds  and  currents  used  to  drive  the 
models  during  the  time  periods  of  the  AIDJEX  experiment  since  this  would 
require  exact  knowledge  of  the  winds  and  currents  everywhere.  It  is 
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interesting,  however,  to  evaluate  the  wind  and  current  errors  necessary 
to  produce  the  observed  errors. 

McPhee  (1980)  has  shown  that  a  free-drift  model  simulates  ice  drift 
accurately  in  the  summer  (July  through  mid-September).  Errors  in  ice 
velocities  may  thus  reasonably  be  ascribed  to  wind  and  current  errors 
for  this  time.  In  Figure  5. 13  mean  velocity  errors  and  standard  deviations 
of  errors  for  free-drift  are  plotted  as  a  function  of  time  of  year.  In 
the  summer  the  mean  daily  velocity  error  is  roughly  0.01  m/s  (~1  km/day) 
with  a  standard  deviation  of  0.03  m/s  (-3  km/day).  For  long-term  free- 
drift,  the  ice  velocity  relative  to  geostrophic  currents  is  essentially 
a  linear  function  of  the  geostrophic  wind  speed  U  : 

v  -  v  «  N  **  U  (5.29) 

g  a 

h  h 

where  N2  -  0.014.  Thus,  if  U  and  N  are  perfectly  accurate,  errors 
a  a 

in  ice  velocity  will  equal  errors  in  geostrophic  current  determinations. 

The  geostrophic  currents  in  the  Beaufort  Sea  range  from  0.01  to  0.03  m/s 
and  are  not  well  known  (Hunkins,  1974).  Highly  energetic  eddy  currents 
with  speeds  up  to  0.60  m/s  and  diameters  of  10  to  20  km  may  further 
increase  the  variability  in  currents  (Hunkins,  1980).  Errors  in  geostrophic 
currents  can  thus  entirely  account  for  the  errors  and  standard  deviation 
of  errors  in  ice  velocity.  Improvement  of  these  errors  is  unlikely 
unless  some  better  measure  of  dynamic  topography  of  the  Beaufort  gyre 
becomes  available. 

Suppose  instead  that  the  geostrophic  currents  are  perfectly  known 

l< 

and  that  the  wind  is  the  only  source  of  error.  If  N  is  constant, 

a 

then  Equation  (5.29)  implies  errors  of  0.7  m/s  in  the  geostrophic  wind 
with  a  standard  deviation  of  2.1  m/s.  Wind  errors  of  about  2  m/s  are 
not  unusual;  the  mean  geostrophic  wind  in  summer  is  ~10  m/s  (Albright, 

1980),  so  this  represents  an  error  of  only  20  percent.  Independent 
measurements  of  geostrophic  winds  differ  by  about  13  percent.  Such  an 
error  can  account  for  a  large  part  of  the  errors  in  ice  velocity. 
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Figure  5.13.  Mean  Error  Magnitudes  and  Standard  Deviations  as  a 
Function  of  Time  of  Year 
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Denner  and  Ashim  (1980)  have  analyzed  errors  between  AIDJEX  winds 
and  those  produced  by  Fleet  Numerical  Weather  Central  (FNWC) .  For  winds 
greater  than  5  m/s  the  mean  error  is  ~40  percent.  Variation  of  this 
_e rror  is  not  given.  Use  of  FNWC  pressure  maps  to  predict  ice  trajectories 
should  therefore  give  rise  to  errors  4  times  as  great  as  those  produced 
using  AIDJEX  pressure  data. 

The  standard  deviations  calculated  here  may  be  used  to  predict 

errors  in  future  free-drift  calculations.  As  seen  in  Figure  5.13,  the 

error  magnitudes  are  uniformly  much  smaller  than  the  standard  deviations. 

This  leads  to  the  possibility  that  these  errors  are  not  signif icantly_ 

different  from  zero.  The  relevant  statistic  to  check  is  t  ,  =  ‘ ^ 

„  n—  1  S 

where  n  =  number  of  days  in  comparison,  |e|  =  average  error  magnitude 
and  S  =  standard  deviation.  If  the  error  is  95  percent  significant, 
t  ^  t  wHpitp  t? 

n-l  0.95,n-l  0.95,n-l  is  the  students  t-distribution 

(Morrison,  1976,  p.  367).  For  most  of  the  calculations  using  the 
free-drift  model,  the  error  is  not  significant.  For  75RUN1F  and  75RUN2B, 
a  significant  error  is  apparent,  but  only  for  75RUN2B  is  it  unequivocal; 

5.3.2  Seasonal  Variations  in  Performance 


Ice  conditions  in  the  Arctic  change  considerably  with  season.  We 
expect  to  see  this  reflected  in  the  performance  of  the  two  models 
presented.  Figure  5.13  shows  the  magnitudes  and  standard  deviations  of 
velocity  errors  described  in  the  previous  sections  as  a  function  of  time 
of  year.  The  error  magnitudes  of  free-drift  calculations  appear  to  show 
a  minimum  steady  value  of  from  0.006  to  0.010  m/s  from  May  through 
December  with  larger  errors  for  February  and  April.  This  might  be 
deceptive.  The  three  calculations  for  February,  April  and  May  are  made 
from  data  sets  of  10  to  15  days,  while  the  calculations  for  June  through 
December  include  30  days  worth  of  data  each.  The  central  limit  theorem 
holds  that  the  observed  mean  error  should  approach  the  true  error  magnitude 
as  the  number  of  observations  increases.  Calculations  based  on  short 
time  periods  may  thus  be  expected  to  estimate  mean  errors  less  accurately. 
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Nonetheless,  we  do  know  that  the  free-drift  model  should  work  better  in 
the  summer  since  the  ice  pack  has  essentially  no  strength  at  this  time. 

The  large  errors  in  February  and  April  reflect  this  fact. 

Standard  deviations  of  velocity  errors  show  a  much  more  pronounced 
seasonal  effect  which  appears  to  be  independent  of  sample  size.  The 
lowest  values  occur  from  May  through  September  ranging  from  0.019  to 
0.035  m/s.  The  standard  deviation  increases  in  the  fall  to  a  maximum  of 
0.095  m/s  in  February  and  drops  again  in  the  spring.  Mean  wind  speeds 
in  the  Arctic  do  not  show  a  similar  variation  from  season  to  season, 
so  the  increased  variation  in  the  winter  is  not  due  to  increased  wind 
input.  Free-drift  simply  does  not  work  well  in  the  winter  because  the 
ice  strength  becomes  appreciable  at  this  time.  When  the  ice  strength  is 
accounted  for,  by  the  complete  ice  dynamics  model,  the  standard  deviation 
in  February  drops  to  0.036  m/s,  which  is  comparable  to  the  summer  free- 
drift  value. 

5.3.3  Accumulating  Errors  in  Trajectories 

Over  a  period  of  days  or  weeks  it  is  important  to  know  how  trajectory 
errors  accumulate  and  how  much  variation  can  be  expected  in  the  errors. 

The  trajectory  error  vector  simply  grows  linearly  with  the  number  of 
days: 

*  Ne  .  (5.30) 

The  variance  also  grows  linearly  as  in  a  random  walk: 

=*  NS*  (5.31) 

2 

where  S£  is  the  variance  of  daily  velocity  errors.  The  standard 
deviation  thus  grows  as  the  square  root  of  N 


(5.32) 
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In  practice,  errors  can  accumulate  in  two  ways,  depending  on  whether 
the  direction  of  the  mean  error  is  known  reliably  in  addition  to  its 
magnitude.  Figure  5.14  shows  a  schematic  of  a  modeled  ice  trajectory. 

In  .this  case,  the  mean  daily  error  is  assumed  to  be  a  known  vector.  The 
length  of  the  predicted  error  vector  grows  linearly  with  the  number  of 
days  in  the  model.  Each  circle  represents  one  standard  deviation  in 
error  magnitude;  it  encloses  the  area  within  which  the  probability  of 
finding  the  true  ice  position  is  0.683.  The  error  is  assumed  to  have  a 
normal  Gaussian  distribution  about  its  mean  so  radii  of  different  pro¬ 
babilities  may  be  easily  found  by  solving 


for  r  where  p(r) 


is  the  desired  probability. 


(5.33) 


In  the  context  of  searching  for  an  oil  spill  whose  trajectory  has 
been  predicted  by  the  free-drift  model,  note  that  the  area  to  be  searched 
with  a  given  probability  of  success  increases  linearly  with  time.  If 
the  direction  of  the  error  vector  is  unknown,  the  area  to  be  searched 
increases  substantially.  The  search  area  corresponds  to  an  envelope 
about  a  set  of  circles  of  radius  whose  centers  lie  on  a  circle  of 

radius  |£t|  ,  Figure  3.15.  These  describe  a  circle  of  radius 
r  *  l£tl  +  St  to  be  searched.  Since  | e 1 1  increases  linearly  with 
time,  the  search  area  from  this  component  of  t  increases  as  the  square 
of  time.  This  is  a  much  more  realistic  situation  since  we  do  not  have  a 
good  measure  of  mean  error  vectors  averaged  over  long  time  periods.  If 
these  were  available  they  could  be  used  to  correct  the  model,  probably 
by  correcting  the  geostrophic  current,  and  the  model  would  have  a  zero 
mean  error.  From  Figure  5.13  it  is  apparent  that  a  mean  free-drift 
error  magnitude  of  -0.01  m/s  is  typical  from  May  through  December  1975. 
The  standard  deviation  for  this  time  period  is  about  0.03  m/s.  After  7 
days,  this  corresponds  to  a  search  radius  of  13  km.  After  1  month, 
the  radius  would  be  46  km.  This  compares  very  favorably  with  the  25- 
year  variability  of  motions  which  could  locate  an  oil  spill  to  within  an 
ellipse  50  km  wide  by  150  km  long.  The  area  to  be  searched  is  thus 
about  one-third  that  required  if  no  trajectory  is  calculated. 
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Figure  5.14.  Accumulation  of  Errors  in  Trajectory  when  Mean  Error  Vector  Is 
Known.  Circles  Represent  the  Area  to  Be  Searched  with  a 
Probability  of  0.683  of  Finding  the  True  Trajectory  After  N  Days. 
Radii  Are  Given  by  St  =s/N  S£.  Error  Vectors  are  Accumulated 
Mean  Errors  Given  by  et  =  Nf 
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Figure  5.15.  Accumulation  of  Errors  in  Trajectory  if  Only  Error  Magnitude  and 
Standard  Deviation  Are  Known.  Circles  Represent  the  Area  to  Be 
Searched  with  a  Probability  of  0.683  of  Finding  the  True  Trajectory 
After  N  Days.  Radii  Are  Given  by  (tt|  +  St 
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From  October  to  December  1975,  the  free-drift  standard  deviation 
doubles,  but  this  does  not  cause  a  significant  increase  in  the  search 
radius  since  the  mean  error  remains  low.  Assuming  S£  =  0.07  m/s  and 
| e[  =  0.01  m/s  gives  a  search  radius  of  68  km. 

In  the  winter,  the  free-drift  mean  error  can  triple.  Assuming  a 
mean  error  of  0.03  m/s  and  a  standard  deviation  of  0.07  m/s  gives  a 
search  radius  of  130  km  which  is  comparable  to  the  natural  variability 
of  motion.  During  this  time,  a  free-drift  calculation  performs  only  a 
little  better  than  assuming  that  the  ice  does  not  move  at  all. 

The  complete  ice  dynamics  model,  on  the  other  hand,  does  quite  well 
in  the  winter.  If  a  mean  error  of  0.005  m/s  and  a  standard  deviation  of 
0.036  m/s  are  representative,  this  gives  a  search  radius  of  only  35  km. 


Chapter  6 

RESOURCE  REQUIREMENTS 


6.1.  GENERAL  DISCUSSION 

All  calculations  performed  to  date  by  the  authors  have  used  a 
CDC  6400  computer.  Resource  requirements  reflect  this  fact.  If  another 
computer  is  used,  then  some  of  the  requirements,  central  processing  time 
in  particular,  can  change  dramatically .  The  user  can  expect  comparative 
increases  or  decreases  for  both  models.  It  should  also  be  noted  that 
computation  times  and  other  resource  requirements  could  be  decreased 
substantially  by  modifying  the  codes.  This  is  especially  true  of  the 
complete  ice  dynamics  model  which  has  been  coded  to  allow  great  flexibility 
because  it  is  a  research  tool,  not  an  operational  forecasting  system. 

Estimates  of  resource  requirements  are  given  separately  for  each 
model.  These  estimates  can  be  affected  strongly  by  the  format  of  input 
data  and  by  the  format  of  the  output  and  the  printer/plotter  device  used 
to  present  it.  We  assume  that  an  operational  forecasting  system  will 
receive  wind  predictions  from  a  U.S.  governmental  agency,  such  as  the 
NOAA  National  Weather  Service  or  the  Navy  Fleet  Numerical  Oceanographic 
Center,  or  the  Canadian  Atmospheric  Environment  Service  in  a  form  that 
is  immediately  useful.  A  similar  assumption  is  made  for  any  ocean 
current  data  that  are  available.  Also,  we  assume  that  a  line  printer 
and  a  plotter  (either  an  x-y  plotter  or  a  CRT  display  with  hard  copies) 
are  available.  Plotting  software  is  also  required  to  allow  a  concise 
display  of  ice  motion  fields  and  trajectories.  These  requirements  are 
similar  for  both  the  free-drift  and  complete  ice  dynamics  models. 

6.2  FREE-DRIFT  MODEL 


6.2.1  Computer 


The  free-drift  computer  program  determines  the  ice  velocity  at  a 
location  independently  of  surrounding  conditions.  Therefore,  except  for 
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the  need  to  have  wind  and  current  fields  available,  the  trajectory  of  a 
given  ice  floe  (a  material  point  in  the  mathematical  model)  is  found 
independently  of  surrounding  locations.  This  fact  implies  that  computer 
memory  and  timing  requirements  are  minimal. 

The  computer  program  requires  about  40  K  octal  words  during  compilation. 
Execution  is  estimated  to  require  only  a  few  thousand  words,  say  3  K  octal. 
This  does  not  include  storage  arrays  for  winds  and  currents.  There  are 
no  requirements  for  mass  storage  units  such  as  disks  or  tape  drives. 

Central  processing  unit  timing  requirements  are  roughly  3  ms  per 
velocity  estimate.  There  is  no  time  step  limitation,  so  free-drift 
velocity  estimates  can  be  made  each  6  hours  or  whenever  wind  and  current 
data  are  available.  These  execution  costs  are  small  enough  to  be  considered 
negligible.  Similarly,  the  peripheral  processor  and  input/output  times 
are  negligible. 


It  is  expected  that  the  free-drift  model  could  be  executed  on  the 
smallest  possible  minicomputer,  and  probably  on  a  hand-held  programmable 
calculator.  No  attempts  have  been  made  to  do  so,  however.  It  is  important 
to  be  able  to  display  ice  trajectories  on  a  CRT  or  x-y  plotter  to  give 
confidence  to  the  user.  These  plotting  requirements  overshadow  actual 
needs  for  computing  ice  motions. 


6.2.2  Manpower 


An  estimate  of  ice  motion  during  one  time  increment  requires  initial 
position,  average  wind  and  current  data  and  ice  thickness.  The  manpower 
required  to  input  these  data  and  execute  the  computer  program  is  too 
small  to  estimate  accurately,  whether  the  simulation  is  being  performed 
for  the  first  time  or  during  an  ongoing  operation.  Efforts  to  produce 
wind  and  current  input  have  been  considered  separately,  as  have  the 
plotting  requirements.  These  manpower  requirements  will  be  dominated  by 
the  need  to  have  personnel  available  full-time  to  calculate  and  interpret 
results . 


6.3  COMPLETE  ICE  DYNAMICS  MODEL 


6.3.1  Computer 

All  simulations  performed  to  date  by  the  authors  have  used  numerical 
grids  roughly  12  x  15  cells  in  each  direction,  about  208  nodes.  Such  a 
grid  can  cover  the  Beaufort  Sea  with  cells  of  about  A  =  40  km  .  If  the 

same  region  is  covered  with  more  cells  of  smaller  sizes,  then  the  computing 

_3 

time  is  proportional  to  A  .  This  occurs  because  the  number  of  nodes 

-2 

is  proportional  to  A  and  the  maximum  time  step  is  proportional  to 
A  .  Thus,  more  cycles  are  needed  to  advance  the  solution  the  same 
increment  in  time. 

The  computer  program  requires  approximately  100  K  octal  words 
during  execution.  This  includes  storage  for  both  program,  data  and 
output  arrays.  Central  memory  is  not  increased  appreciably  if  more 
cells  are  used  because  solutions  are  stored  on  disks  or  other  mass 
storage  units  during  computations.  A  disk  (or  other  mass  storage  unit) 
is  required  to  provide  intermediate  storage  of  solutions.  Two  time 
levels  of  solutions  are  stored  simultaneously.  Each  solution  describes 
about  26  variables  for  each  mode.  This  is  affected  by  the  number  of 
thickness  levels  used  to  describe  the  ice  thickness  distribution.  A 
tape  drive  unit  is  used  only  as  a  backup  for  permanent  file  storage. 

Central  processing  unit  (CPU)  timing  requirements  are  roughly 
20  ms  per  cell  per  cycle,  exclusive  of  printout  and  plotting  operations. 
This  CPU  time  breaks  down  into  about  14  percent  for  constitutive  law, 

35  percent  for  thickness  distribution,  31  percent  for  momentum  balance 
and  the  remaining  20  percent  for  array  shifting  and  other  overhead 
functions.  Time  steps  are  on  the  order  of  3  minutes,  so  approximate  ly 
480  cycles  are  needed  to  increment  the  solution  forward  1  day.  This  is 
the  smallest  time  step  used  to  date  in  all  simulations.  Peripheral 
processing  times  are  observed  to  be  about  one-quarter  to  one-half  of  CPU 
times.  Input/output  times  are  observed  to  be  about  one-third  to  two-thirds 
of  CPU  times.  This  last  estimate  depends  linearly  on  the  frequency  at 
which  solutions  are  stored  for  later  printing  and  plotting. 
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The  computer  resources  required  to  execute  the  complete  ice  dynamics 
model  can  vary  by  a  factor  of  4  to  8,  either  longer  or  shorter  depending 
on  many  specific  parameters.  A  users'  guide  prepared  by  D.  R.  Thomas 
(unpublished  manuscript)  has  indicated  how  such  choices  can  affect 
timing.  As  an  example  of  actual  costs,  the  17-day  simulation  of  winter 
ice  dynamics  presented  by  Pritchard  (1981)  required  approximately  22,600 
CPU  seconds  to  execute.  The  grid  included  12  x  16  =  192  cells  and  a 
time  step  of  180  seconds  per  cycle  was  used.  This  time  is  about  three- 
fourths  of  the  average  estimate  presented  in  this  section.  For  the  sake 
of  completeness,  it  is  noted  that  known  modifications  could  reduce  this 
cost  by  an  estimated  factor  of  about  8,  or  to  12  percent  of  the  stated 
cost . 

6.3.2  Manpower 


An  individual  simulation  of  sea  ice  dynamics  using  the  complete 
model  requires  manpower  to  set  up  and  initialize,  to  execute  and  to 
print  or  plot  output  and  interpret  results.  The  problem  set  up  and 
initialization  requires  that  a  grid  be  selected  and  positions  of  all 
nodes  be  specified.  Initial  ice  velocity,  stress  and  thickness  distribution 
are  required  throughout  the  region.  Only  the  thickness  distribution 
needs  to  be  input  with  care  at  the  beginning  of  a  new  simulation  (or 
forecast).  The  grid,  velocity  and  stress  can  be  input  once  for  any 
region  of  interest.  Execution  of  the  computer  program  typically  is 
accomplished  by  a  computer  facility  staff  available  as  part  of  the 
computing  cost.  Plotting  of  output  data  and  interpretation  of  results 
depend  on  the  purpose  of  a  simulation.  If  a  basic  understanding  of  ice 
dynamics  in  a  new  region  or  at  a  different  time  of  year  is  needed,  then 
interpretation  of  a  simple  month-long  simulation  is  a  two-to-four  man- 
month  effort.  As  a  standard  forecast,  this  postprocessing  is  minimal. 

For  an  operational  forecasting  system,  manpower  estimates  are  expected 
to  be  roughly  comparable  to  free-drift  model  operations. 
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Chapter  7 

CONCLUSION  AND  RECOMMENDATIONS 

„  The  free-drift  model  has  been  tested  as  a  function  of  season  and 
geographic  location.  Errors  in  free-drift  calculations  have  been  compared 
with  the  natural  variability  of  ice  motion  and  with  the  more  sophisticated 
complete  ice  dynamics  model. 

The  most  significant  error  is  the  accumulated  trajectory  error. 

This  determines  the  area  to  be  searched  (about  the  predicted  position) 
after  a  given  number  of  days  for  the  true  location  of  a  point  (oil 
spill)  on  the  ice  pack.  After  a  month,  the  ice  can  move  considerable 
distances  (50  to  150  km)  in  the  summer  or  winter.  The  free-drift  model 
simulates  a  1  month  trajectory  to  within  50  km  in  the  summer.  Considering 
the  low  cost  and  ease  of  implementing  the  free-drift  model,  this  represents 
a  significant  improvement  over  a  search  from  the  point  of  origin.  In  the 
winter,  free-drift  appears  to  offer,  at  best,  only  a  marginal  improvement 
over  searching  the  historical  variability  ellipse.  The  complete  ice 
dynamics  model  does  much  better  in  the  winter  and  is  comparable  in 
accuracy  to  free-drift  in  the  summer.  This  model  is  more  difficult  to 
implement,  however,  requiring  tracking  buoys  on  the  perimeter  of  the 
area  and  accurate  winds  over  the  entire  area.  The  free-drift  model,  in 
comparison,  requires  only  local  winds  near  the  point  of  interest.  The 
additional  costs  incurred  in  implementing  the  complete  model  must  be 
weighed  against  the  improvement  in  predicted  motions. 

Both  the  free-drift  and  the  complete  ice  dynamics  models  would 
benefit  substantially  from  an  improved  geostrophic  current  field  in  the 
Beaufort  Sea.  Such  a  current  field,  if  accurately  determined,  could 
reduce  mean  daily  errors  to  zero,  though  variability  in  errors  would 
remain.  Free-drift  trajectory  errors  in  the  summer  should  fall  to  less 
than  20  km  and  in  the  winter,  to  less  than  40  km.  The  free-drift  model 
could  easily  be  used  to  obtain  the  geostrophic  current  for  any  summer  in 
which  tracking  buoy  data  are  available.  The  ice  velocity  field  given 
by  the  free-drift  model  is  velocity  relative  to  geostrophic  current. 
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Monthly  mean  differences  between  modeled  motions  and  observed  motions 
thus  give  the  geostrophic  current  directly.  These  currents  should  vary 
slowly  over  periods  of  months  or  longer  and  could  also  be  used  to 
improve  the  accuracy  of  the  complete  ice  dynamics  model  in  the  winter. 
Accurate  currents  should  reduce  monthly  complete  model  trajectory  errors 
in  the  winter  to  less  than  20  km. 

We  recommend  that  the  free-drift  model  be  implemented  to  aid  in 
tracking  oil  spills  from  May  until  October  and  a  complete  ice  dynamics 
model  be  used  during  the  winter  and  early  spring  months.  The  free- 
drift  model  should  also  be  used  to  construct  a  geostrophic  current  field 
for  the  Beaufort  Sea. 
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